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ABSTRACT 


This  report  describes  23  investigations  in  the  field  of  seismic 
discrimination.  These  are  grouped  as  follows:  estimation  of 

magnitude  by  networks  and  single  stations  (6  contributions); 
studies  relating  to  earth  structure  and  scattering  processes 
(7  contributions);  studies  in  focal  depth  determination,  source 
mechanism  discrimination,  and  general  seismology  (6  contri- 
butions); and  recent  developments  in  computer  systems  and 
software  for  seismic  data  processing  (4  contributions). 
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SUMMARY 


This  is  the  twenty -fourth  Semiannual  Technical  Summary  report  describing  the  activities 
of  Lincoln  Laboratory,  M.I.T.,  in  the  field  of  seismic  discrimination.  These  activities  involve 
research  into  the  fundamental  seismological  problems  associated  with  the  detection,  location, 
and  identification  of  earthquakes  and  nuclear  explosions.  We  also  are  concerned  with  the  devel- 
opment of  methods  for  the  handling  and  analysis  of  large  quantities  of  global  seismic  data,  and 
the  application  of  these  methods  to  data  management  system  design  and  the  optimum  extraction 
of  scientific  information  from  high-quality  digital  data. 

A detailed  investigation  into  the  problems  associated  with  the  estimation  of  magnitude  by  a 
network  of  seismic  stations  has  been  initiated.  This  was  prompted  by  some  earlier  studies  that 
suggested  that  the  effects  of  station  detection  thresholds  may  be  to  insert  a positive  bias  into 
network  magnitudes  for  small  events.  The  problem  is  primarily  a statistical  one,  and  a series 
of  statistical  models  have  been  developed  in  order  to  facilitate  the  removal  of  these  effects. 

The  principal  features  of  these  models  are  described  in  this  report.  The  next  important  problem 
is  the  optimum  determination  of  the  parameters  describing  the  detection  characteristics  and 
bias  of  each  station  in  a network.  Some  preliminary  results  along  these  lines  are  presented  in 
this  report.  Other  studies  in  magnitude  estimation  have  concentrated  on  existing  bulletin  data, 
and  several  investigations  show  that  current  methods  of  determining  magnitude,  particularly 
body  wave  magnitude  m^,  are  open  to  some  criticism. 

We  have  continued  to  attack  the  general  problem  of  earth  structure  and  heterogeneity.  An 
inversion  of  free  oscillation  data  to  obtain  the  attenuation  parameter  Q shows  rather  low  values 
of  Q in  both  upper  and  lower  mantles,  and  no  evidence  for  a frequency  dependence  of  Q.  Earlier 
work  on  tracing  the  paths  of  Rayleigh  waves  across  the  Eurasian  continent  has  been  extended  by 
the  inclusion  of  more  information  on  crustal  structure.  Severe  deviations  from  great-circle 
paths  and  large  amplitude  anomalies  are  predicted,  and  the  agreement  with  observation  is  good. 
Studies  of  the  P-wave  codes  from  Novaya  Zemlya  explosions,  using  wave  number  analysis  and 
adaptive  deconvolution,  show  complex  arrivals  generated  in  the  vicinity  of  the  source.  Other 
studies  described  include  reflections  from  the  600-km  discontinuity,  and  the  applications  of 
spatial  correlation  and  transmission  holography  to  seismic  data. 

Section  III  of  this  report  describes  a series  of  investigations  related  to  seismic  discrimina- 
tion and  general  seismology.  Further  investigations  of  the  Maximum  Entropy  method  of  spectral 
analysis  have  focused  on  the  determination  of  the  order  of  a time  series,  and  the  optimization  of 
the  length  of  the  prediction  error  filter.  Application  of  the  method  to  the  network  determination 
of  focal  depth  is  included.  Effects  of  source  geometry  on  the  discrimination  between  events  with 
different  source  mechanisms  are  being  studied.  A series  of  travel-time  tables  based  on  a para- 
metric representation  are  described.  These  have  the  advantages  of  high  accuracy  and  low 
computer -storage  requirements.  An  analysis  is  described  which  links  the  concept  of  the  moment 
tensor  of  an  earthquake  with  its  radiated  energy.  An  application  of  normal  mode  theory  leads 
to  a new  method  for  the  calculation  of  the  permanent  displacement  fields  of  an  earthquake. 

•We  have  continued  to  develop  the  software  and  systems  necessary  for  the  analysis  of  digital 
seismic  data  from  a global  network.  We  have  adapted  previous  software  into  a form  suitable 


for  the  display  and  processing  of  SRO  data.  An  interactive  system,  utilizing  the  IBM  370/168 
at  Lincoln  Laboratory,  has  been  designed  for  the  comparison  of  observed  waveforms  with  those 
predicted  by  source  mechanism  theory.  We  are  continuing  to  update  the  seismic  documentation 
files  that  we  constructed  using  the  NLS  system.  Also,  we  are  developing  the  additional  software 
necessary  for  the  full  utilization  of  our  connection  to  the  ARPANET  and  the  datacomputer  at  the 
Computer  Corporation  of  America. 


M.  A.  Chinnery 


GLOSSARY 


ARPANET 

DARPA  Computer  Network 

CCP 

Control  and  Communication  Processor 

DARPA 

Defense  Advanced  Research  Projects  Agency 

EDR 

Earthquake  Data  Reports 

HGLP 

High  Gain  Long  Period  (Network) 

ISC 

IWWSS 

International  Seismological  Center 
Integrated  World-Wide  Seismic  System 

LASA 

Large  Aperture  Seismic  Array 

NEP 

NIC 

NLS 

NORSAR 

NTS 

Network  Event  Processor 
Network  Information  Center 
On-Line  System 
Norwegian  Seismic  Array 
Nevada  Test  Site 

PDE 

PEM 

Preliminary  Determination  of  Epicenters 
Parametric  Earth  Model 

SATS 

SIP 

S/N 

SRO 

sus 

Semiannual  Technical  Summary 
Seismic  Information  Processor 
Signal  to  Noise  (Ratio) 

Seismic  Research  Observatory 
Seismic  User  Subsystem 

USCGS 

USGS 

U.S.  Coast  and  Geodetic  Survey 
U.S.  Geological  Survey 

VESPA 

Velocity  Spectral  Analysis 

WWSSN 

World-Wide  Standard  Seismograph  Network 
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SEISMIC  DISCRIMINATION 


I.  ESTIMATION  OF  MAGNITUDE 

A.  PROBLEMS  IN  MAGNITUDE  ESTIMATION 

Lincoln  Laboratory  has  embarked  on  an  in-depth  investigation  of  the  problems  involved  in 
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the  estimation  of  magnitude,  using  both  single  stations  and  networks.  Previous  Lincoln  studies, 
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together  with  more  recent  work  of  Ringdal,  have  demonstrated  the  importance  of  including  sta- 
tion detection  thresholds  and  biases  when  attempting  to  estimate  the  magnitude  of  small  events 
using  a network.  We  have  extended  these  results,  and  have  begun  to  apply  them  to  network  data. 

The  basic  problems  can  be  divided  into  two  separate  cases.  The  first  problem  concerns 
the  development  of  methods  for  the  estimation  of  magnitude  by  an  idealized  network.  The  cru- 
cial point  here  is  that  a report  from  each  station  in  the  network  is  received  for  each  event.  This 
report  will  consist  of  either  "seen,"  when  an  amplitude  and  period  will  be  recorded,  or  "not- 
seen,"  when  the  instrument  noise  level  will  be  recorded.  Such  an  idealized  network  is  not  pres- 
ently available,  but  the  new  ARPA  global  network  of  arrays,  SROs,  and  HGLP  stations  may  be 
used  in  this  way. 

The  second  class  of  problem  involves  the  data  from  imperfect  networks,  such  as  the  group 
of  stations  currently  used  to  produce  the  ISC  and  PDE  Bulletins.  In  these  cases,  the  "not  seen" 
information  is  generally  unreliable,  due  to  inconsistencies  in  both  instrumental  status  and  oper- 
ator characteristics. 

In  Secs.  B through  F that  follow,  we  outline  a series  of  statistical  models  that  may  be  appli- 
cable to  the  above  problems,  and  explore  ways  of  estimating  station  detection  threshold  and  bias. 
Some  preliminary  results  from  the  application  of  these  models  are  included. 

Perhaps  inevitably,  detailed  examination  of  these  problems  has  raised  a series  of  new  ques- 
tions that  were  not  anticipated.  Aspects  of  this  work  currently  under  investigation,  and  which 
will  be  reported  in  later  publications,  are  as  follows: 

(1)  Presumably  clipping,  and  the  difficulty  in  measuring  very  large  swings 
on  a photographic  seismogram,  leads  to  effects  very  similar  to  a detec- 
tion threshold.  Can  they  be  described  by  similar  statistical  models? 

And  how  does  this  affect  the  magnitude  estimation  of  large  events  in 
existing  catalogues? 
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(2)  To  date,  following  Ringdal,  we  have  assumed  the  frequency-magnitude 

relation  to  be  linear.  In  the  case  of  mK,  this  is  a questionable  assump- 
3 D 

tion.  Newer  models  containing  a quadratic  relation  are  being  developed 
and  applied. 

(3)  One  contribution  to  the  variation  in  station  amplitude  A observed  over  a 
network  for  a given  event  is  due  to  path  differences,  or  Q.  However,  Q 
will  also  affect  the  observed  dominant  period  T.  Recording  log  A/T  only 
partially  corrects  for  this  effect  (see  Sec.  F).  Also,  this  raises  the  com- 
plication that  instrument  response  should  be  removed  from  station  read- 
ings before  the  station  magnitude  is  computed,  and  this  is  not  a routine 
correction  at  all  stations. 
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It  appears  likely  that  it  will  never  be  possible  to  completely  remove  these  type  of  difficul- 
ties from  existing  earthquake  catalogues.  However,  it  may  be  possible  to  remove  systematic 
biases  from  these  data  sources,  and  thereby  reduce  the  observed  scatter  in  observations. 

M.A.  Chinnery 


B.  STATISTICAL  MODELS  FOR  MAGNITUDE  ESTIMATION 

4 12 

We  have  extended  the  earlier  work  of  Freedman,  Kelly  and  Lacoss,  and  Ringdal,  in  order 
to  account  for  the  detection  characteristics  of  seismic  stations  and  networks.  We  first  list  the 
assumptions  and  definitions  used  in  the  models,  and  then  outline  the  models  themselves  for  a 
number  of  useful  cases. 


Assumptions  and  Definitions 

1.  When  the  i ^ station  in  a network  observes  an  event  with  true  magnitude  m,  it  will  re- 
cord a value  of  [log(A/T)J,  where  A is  the  amplitude  and  T the  period  of  the  signal.  For  con- 
venience, we  set 

y.=[log|].  . (1-1) 

We  define  the  station  magnitude  m.  as 

m.  = y.  + Qi  (1-2) 

where  Q.  is  the  station  distance-depth  correction  for  the  event,  at  the  i**1  station.  We  then  as- 
sume that  m.  is  related  to  m by  an  expression  of  the  form 

mi  = m + B.  + €.  (1-3) 

where  B.  is  the  station  bias,  and  e.  is  a Gaussian  random  variable  with  zero  mean  and  variance 
2 1 1 
ai  ’ 

2.  We  assume  that  the  probability  that  an  event  will  be  detected  at  the  station  is  given  by 

Pr  {detection/y.}  = f 1 - — exp  {-1/2  [(z  - G.)/y.]2 3}  dz 

1 J— ° o/2?  y.  1 


$ 


(1-4) 


Gj  is  the  detection  threshold  for  observations  of  y.,  and  corresponds  to  a 50-percent  detection 
probability,  y.  determines  the  shape  of  the  detection  probability  curve. 

3.  In  certain  cases,  it  will  be  more  convenient  to  express  the  detection  characteristics  in 
terms  of  the  station  magnitude  iru.  We  may  then  attempt  to  represent  the  detection  probability 
curve  by  an  expression  of  the  form 


/m.  — G.  \ 
Pr  {detection/m.}  = — -J 

\ 1 i ' 


(1-5) 


Equation  (1-5)  is  formally  equivalent  to  Eq.  (1-4),  with  y?  = y^  and  G?  = G.  + Q.,  if  all  the  events 
considered  are  located  within  a small  region,  so  that  paths  and  Q values  to  the  network  are  the 
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same  for  all  events.  If  the  events  lie  within  a distributed  source  area,  Eq.  (1-5)  is  an  approxi- 
mation, and  G?  and  y?  will  in  general  have  no  simple  relation  to  G.  and  y 

4.  We  assume  that  the  true  magnitudes  of  earthquakes  obey  the  Gutenberg- Richter  rela- 
tionship, in  the  form 

InN  = a —0m  (1-6) 

where  N is  the  number  of  events  in  a small  magnitude  interval  centered  at  m,  and  a and  p are 
constants.  It  is  necessary  that  p have  the  same  value  at  all  parts  of  the  source  region,  for  if 
this  is  not  so,  or  if  the  magnitude -frequency  relation  is  not  linear,  new  complications  arise. 
These  will  be  the  subject  of  future  studies. 

5.  In  the  presence  of  a detection  threshold,  a station  may  report  either  a value  of  y.  (or 
m^),  or  the  report  "not  seen."  In  order  to  include  a proper  statistical  treatment  of  the  last  pos- 
sibility, it  is  necessary  to  introduce  new  variables  ai  and  \x . defined  as  follows: 


a. 

l 


"not  seen" 


if  event  detected 
if  event  undetected 


mi 


"not  seen" 


if  event  detected 
if  event  undetected 


Statistical  Relationships 

1.  From  assumption  1 above,  the  conditional  distributions  of  y^  and  m.  for  events  with 

true  magnitude  m are: 


g°  {y^m}  dy.  = 0 | 


y.  + Q.  — B.  — m 

Ji  l l 


dy; 


(1-7) 


n im.-B.-mi 

h .{m.lm}  dm.  = 0 ! } dm. 

i 1 I I 1 

o(f-)  =-=! — exp(-l/2(x/<r  )1 2]  . 

v V v.  1 


where 


(1-8) 


(1-9) 


2.  In  the  presence  of  a detection  threshold,  these  equations  require  modification.  Now, 
the  conditional  distributions  of  a.,  and  jjl.  (for  both  seen  and  unseen  events)  are 

f /ai“Gi\  (aj+Qj-B.-mi 

— — — J 0 j J da.  (for  "seen”  events) 


g {a.  | m}  da.  = 


-(B.  +m-Q. -G.) 


(1-10) 


(for  "unseen"  events) 
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and 


l^i-Gfl  IM^i  — Bi  — 

| yf  ' j 9 | o” j d^i  (for  "seen"  events) 


h {(ij  | m}  djx.  = 


~(Bi  +m  - G?) 


I F^i 


*2 


(I-H) 


(for  "unseen"  events) 


The  parallelism  of  the  formulations  in  terms  of  y.  and  im  is  clear.  In  the  expressions  that  fol- 
low, we  quote  only  those  for  the  station  magnitudes  m..  Similar  expressions  in  terms  of  y^  are 
easily  derived. 

3.  If  we  consider  only  those  events  "seen"  at  the  i**1  station,  the  conditional  distribution 
for  observed  station  magnitudes  becomes 


H {rm | m}  dm.  = 


$ 


+ m - G* 


dm. 

l 


(1-12) 


The  shape  of  this  distribution  for  the  parameters  m = 3.8,  = 0,  G ? =4.0,  <r.  = 0.3,  and 

yf  - 0.2  is  shown  in  Fig.  1-1.  For  comparison,  the  Gaussian  distribution  with  mean  3.8  and 
standard  deviation  0.3  is  also  shown. 


4.  The  expected  value  of  observed  station  magnitudes  rm,  given  m,  may  be  found  from 
Eq.  (1-12)  to  be 


where 


E {m . | m}  = m + B.  + 


a. 

i 


F 


2 + yf 2 


m + B.  - G.* 


'of  + yf’ 


J ‘ 


(1-13) 


This  equation  is  plotted  in  Fig.  1-2,  with  the  parameters  G?  = 4.0,  B^  = 0,  o\  = 0.3,  and  yf  = 0.2. 

5.  In  general,  if  an  event  of  magnitude  m is  declared  detected  by  a network  of  N stations 
when  it  is  detected  by  at  least  one  station  in  the  network,  the  conditional  distribution  of  the  sta- 
tion observations  p^,  . . . , p^  is 


N 


H {p  | » • • • , M-jsj)  dp  ^ , . . . , dp^j  = 


n h {u.|m}  d,!. 
i=l  1 1 


N 

i - n 

i=l 


- (Ba  + m - Gp 

A2  + y? 


(1-14) 


In  the  case  of  a one-station  "network,"  Eq.  (1-14)  reduces  to  Eq.  (1-12).  In  the  case  of  Eq.  (1-12), 
rm  can  be  used  rather  than  pA  because  the  pi  never  takes  on  the  "not  seen"  value. 
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6.  Introducing  the  frequency-magnitude  relationship,  Eq.(I-6),  we  may  derive  the  distri- 

th 

but  ion  of  y^  and  m . for  those  events  seen  at  the  1 station  without  conditioning  on  m: 


/5(G.-l/2/3y.2-y .)  ( yj  - G. 

Q(yj)  dyj  = /3  e $ 1 


dyi 


(1-15) 


/5(Gf-l/20V*2-m .)  (m.-Gf 

R(rm)  dm.  = 0 e 1 


dm 


i 


(1-16) 


Notice  that  these  distributions  are  independent  of  B.  and  o\. 

7.  Again  using  the  linear  seismicity  assumption,  the  distribution  of  magnitudes  recorded 
by  a network  of  N stations,  given  that  the  event  is  seen  by  at  least  one  station,  becomes 


R(^1, . . . , nN)  dp>1,  . . . , d^N  = 


£«  [ini  h dp.. 


e~^m  dm 


(1-17) 


r 

J—  0 


N 

i - n 
i=l 


-(B.  + m - G.*) 


■FT^ f 


e^m  dm 


8.  The  conditional  distribution  of  the  station  magnitudes  m2  observed  at  station  2,  given 
that  a station  magnitude  has  been  observed  at  station  1,  is 


g {m2|m1}  dm 2 = 


m2  “ B2  + B1  + Pal 


J°l + a 


m2  ~G2 
?2* 


dm. 


4> 


m ^ + B^  B4  G2 


/ 2 , 2" *2 
Ja  1 + <T2  + ?2 


(1-18) 


The  expected  value  of  m2,  given  m^,  is  then 


2 2 
2 ai  + a2 

E {m2|m1}  = mt  + B2  - B1  -jSc^  + 


+<7: 


2 + ^22 


X z 


+ B2  ~B1  ~G2 
> I'  + ff22  + 


(1-19) 


Because  of  the  similar  forms  of  Eqs.  (1-13)  and  (1-19),  Fig.  1-2  may  be  taken  as  a plot  of 
Eq.  (1-19)  if  the  horizontal  axis  is  labeled  m^  and  the  vertical  axis  m2.  The  appropriate  param- 
eters are  now  G^  =4.0,  B2  — B^  — 0a f = 0,  = cr^  = 0.3 /\[T,  and  = 0.2. 


L.  A.  Christoffersson 
R.T.  Lacoss 

M.  A.  Chinnery 
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C.  ESTIMATION  OF  NETWORK  MAGNITUDE  AND  STATION 
DETECTION  PARAMETERS 

Based  on  the  relationships  developed  in  the  previous  section,  a series  of  computer  programs 
have  been  written  for  the  separate  and  joint  estimation  of  network  magnitude  and  the  parameters 
describing  the  detection  characteristics  of  the  stations  within  the  network.  Four  of  these  pro- 
grams, together  with  some  applications,  are  described  below. 


1.  Network  Magnitude  Estimation 

This  program  is  designed  to  make  a maximum-likelihood  estimate  of  the  magnitude  of  an 
event  detected  by  a network,  assuming  that  the  detection  probability  curves  and  biases  are  known 
at  each  station  in  the  network. 

We  assume  that  reports  are  received  from  each  of  the  N stations  in  the  network.  As  in 
the  previous  section,  these  reports  may  consist  either  of  a station  magnitude  m.,  or  the  report 
"not  seen."  The  appropriate  distribution  of  the  is  given  by  Eq.  (1-14),  and  the  log-likelihood 
function  L is  obtained  by  simply  taking  the  logarithm  of  this  equation.  We  obtain 


N 

log  L = E 

i=l 


log  h {ji.  I m}  - log 


N 

i-n 

i=l 


-(B.  + m-G*) 


J- 


«,2 * * * * *  * yf 


(1-20) 


The  program  maximizes  this  function  with  respect  to  m in  order  to  obtain  the  maximum- 
likelihood  estimate  of  m.  B.,  G*,  and  y ? are  assumed  known.  Notice  that,  in  the  case  of  a sin- 
gle event,  we  may  write  G?  = G^  + and  y?  = y.  (see  Sec.B  above). 

The  likelihood  function  (1-20)  differs  slightly  from  that  of  Ringdal,^  since  he  chose  to  include 
in  his  sample  space  the  possibility  that  none  of  the  stations  saw  the  event.  Also,  the  distribution 
h {pj|m}  differs  from  his  since  we  include  the  appropriate  detection  probability. 

There  remain  some  theoretical  difficulties  with  the  approach  outlined  above.  It  is  known 
that  maximum-likelihood  estimates  are  asymptotically  efficient  for  a broad  class  of  distribution 
functions.  However,  we  have  not  yet  been  able  to  demonstrate  that  the  distribution  (1-14)  falls 
into  this  class.  At  present,  this  method  is  the  only  one  that  allows  for  the  inclusion  of  "not  seen" 
information,  and  we  shall  therefore  pursue  these  problems.  Certainly,  initial  results  suggest 
that  the  method  converges  properly  to  reasonable  values,  both  in  this  case  and  the  one  described 
below. 


2.  Joint  Estimation  of  Network  Magnitude,  Station  Biases, 

and  Detection  Thresholds 


Suppose  a set  of  K events  occurs  within  a localized  region,  so  that  the  path  corrections  Q. 

are  common  to  all  events.  We  may  then  form  a total  log-likelihood  function  by  summing  the  log 

likelihoods  for  the  individual  events,  as  given  by  Eq.(I-20).  Noting,  as  before,  that  in  this  case 

Gj*  = G.  + and  y?  = y.,  we  find 


K K 

log  L = Yj  lQg  = E 
j=1  j=1 


N 


E loSh  - log 


i=l 


N 

i - n * 

i=l 


-(Bi  + nij  - G*) 


F * F 


(1-21) 


.th 


where  m.  is  the  true  magnitude  of  the  j event. 
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A program  which  maximizes  this  function  with  respect  to  the  rm,  the  B.,  and  the  has 
been  written,  a.,  y^,  and  Q.  must  be  assumed  known.  It  has  been  tested  on  an  aftershock  se- 
quence of  72  events  located  in  a small  region  south  of  the  island  of  Honshu,  Japan.  Station  mag- 
nitudes were  obtained  from  the  data  base  of  the  International  Seismic  Month5  for  a set  of  15  sta- 
tions that  were  very  carefully  analyzed.  We  have  assumed  that  the  lack  of  a station  magnitude 
corresponds  to  the  verdict  ”not  seen,”  and  this  is  a very  reasonable  assumption  for  this  data  set. 
o\  and  y.  were  preset  to  the  values  0.3  and  0.2,  respectively,  for  all  sfauons.  The  resulting 
estimates  of  station  bias  B^,  and  station  detection  thresholds  and  G?,  are  listed  in  Table  1-1. 
The  extent  to  which  these  estimates  are  affected  by  incorrect  values  of  o.  and  y.  is  not  yet  clear 
(the  inclusion  of  cr.  and  y.  as  free  parameters  in  the  problem  results  in  severe  computational 
difficulties). 

Estimates  of  network  magnitude  may  be  plotted  against  the  usual  definition  of  magnitude 
(the  average  of  the  station  magnitudes  for  those  stations  reporting  the  event).  This  is  shown 
in  Fig.  1-3.  Clearly,  at  low  magnitudes,  the  conventional  estimate  is  biased  high. 


TABLE  1-1 

ESTIMATES^  OF  STATION  BIAS  B., 

1 

AND  STATION  DETECTION  THRESHOLDS  G.  AND 

i 

G * 

i 

Station 

B. 

i 

G* 

i 

Q. 

i 

G. 

i 

LAO 

0.  07 

3.6 

3.7 

-0.  1 

MBC 

0. 29 

4. 1 

3.9 

+ 0.2 

NAO 

0. 00 

3.7 

3.7 

0.0 

RES 

0.38 

4.2 

3.7 

+ 0.5 

HFS 

0.00 

3.8 

3.7 

+ 0. 1 

UBO 

-0.  16 

3.8 

3.7 

+ 0. 1 

KBL 

0.09 

4.  1 

3.9 

+ 0.2 

FFC 

-0.03 

4.3 

3.7 

+ 0.6 

BLC 

0.20 

4.6 

3.7 

+ 0.9 

ALE 

0.05 

4.5 

3.7 

+ 0.8 

FBC 

0.05 

4.5 

3.7 

+ 0.8 

CHG 

-0. 39 

4.0 

3.7 

+ 0.3 

COL 

-0.21 

4.2 

3.8 

+ 0.4 

FCC 

-0. 15 

4.4 

3.7 

+ 0.7 

YKC 

-0.21 

4.4 

3.7 

+ 0.7 

t Results  of  a preliminary  experiment  in  the  application  of  Eq.  (1-21). 
Seventy-two  events  in  an  aftershock  sequence  south  of  Japan  were 
used.  <7  j and  y;  were  assumed  to  be  0.3  and  0.2,  respectively.  Due 
to  the  small  number  of  events  used,  the  standard  errors  of  the  G;  are 
quite  large,  and  values  there  should  be  considered  presently  unreliable. 
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3.  Joint  Estimation  of  Detection  Parameters  and  Seismicity 

Suppose  a record  is  kept  of  all  the  y = log(A/T)  values  recorded  at  a station,  and  the  results 
are  plotted  in  histogram  form.  Assuming  that  the  events  seen  are  drawn  from  a population  dis- 
tributed in  magnitude  according  to  the  Gutenberg- Richter  relation  (1-6),  the  observed  distribu- 
tion of  y values  will  obey  an  equation  of  the  form  of  (1-15).  A log-likelihood  function  may  be 
obtained  from  this  equation: 


log  L = 


£ |log0  + 0(G  - | Py2 

j=l 


-yj)  + log* 


(1-22) 


where  the  summation  is  over  all  the  events  recorded.  A program  has  been  written  to  maximize 
this  likelihood  function  with  respect  to  /3,  G,  and  y.  This  is  equivalent  to  the  method  of  estimat- 
ing /3,  G,  and  y described  by  Kelly  and  Lacoss.1  It  differs  from  the  approach  of  Aki,^  who  as- 
sumed y = 0 (i.e.,  a step  function  detection  probability  curve). 

This  program  has  been  applied  to  228  observations  of  y at  LASA  from  events  in  Japan. 

The  cumulative  distributions  of  observed  and  theoretical  y are  shown  in  Fig.  1-4  which  also 
gives  95-percent  confidence  bands,  based  on  the  Kolmogorov- Smirnov  test.  The  theoretical 
curve  lies  well  within  these  confidence  bands.  The  corresponding  parameter  estimates  are 
0 = 2.12,  G = —0.08,  and  y = 0.19,  with  standard  deviations  of  0.03,  0.02,  and  0.01,  respec- 
tively (p  = 2.12  corresponds  to  a "b-value"  of  0.92).  We  would  expect  these  values  of  G and  y 
to  correspond  to  those  given  in  Table  1-1  for  LAO,  since  both  are  derived  from  events  in  Japan. 
In  fact,  the  agreement  is  better  than  was  expected,  considering  the  difference  in  methods  and 
data. 


4.  Estimation  of  Detection  Parameters  Using  a Pair  of  Stations 

A particularly  promising  method  for  the  estimation  of  detection  parameters  involves  the 
conditional  distribution  of  the  station  magnitudes  m2  at  one  station,  given  the  station  magnitude 
at  a reference  station,  m^.  The  appropriate  distribution  is  given  in  Eq.  (1-18),  and  is  found  to 
be  independent  of  the  detection  characteristics  of  the  reference  station.  A log-likelihood  func- 
tion can  be  constructed  as  above,  and  the  resulting  function  maximized  with  respect  to  the  free 
parameters,  which  in  this  case  are  G^,  and  B = . A program  has  been  written 

to  estimate  these  parameters,  and  is  in  the  final  stages  of  testing. 

Although  a single  station  pair  cannot  permit  the  estimation  of  the  separate  biases  B^  and  B^, 
a method  for  determining  these  over  a network  is  suggested.  First  a reference  station  is  chosen, 
and  then  the  relative  bias  of  all  the  remaining  stations  can  be  found.  If  a new  reference  station 

is  now  chosen  and  the  experiment  repeated,  the  relative  biases  of  all  the  stations  in  the  network 

N 

can  be  found.  Adding  the  usual  condition  that  2 B.  = 0 is  sufficient  to  determine  the  B.  uniquely. 

i=l  1 1 

L.  A.  Christoffersson 
R.  T.  Lacoss 

M.  A.  Chinnery 


D.  STATION  MAGNITUDE  BIAS  AND  ITS  EFFECT 
ON  NETWORK  mb  ESTIMATION 

An  examination  of  individual  station  body  wave  magnitude  (m^)  measurements  has  revealed 
that  these  may  possess  a sizable  bias  relative  to  the  network  m^  (generally  computed  as  an 
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average  of  all  determinations).  The  effect  was  first  noticed  in  a comparison  of  Russian  m,  de- 

7 D 

terminations  with  those  of  the  USCGS  PDE  list,  in  which  it  was  postulated  that  some  stations 
which  report  very  frequently  have  large  negative  biases.  It  is  apparent  that  m^  values  for 
smaller  events,  for  which  only  a few  stations  report  amplitude  measurements,  will  be  more 
greatly  affected  by  such  biases  than  larger  earthquakes,  for  which  many  stations  contribute  to 
the  network  mb  and  thus  average  out  the  effect  of  these  biases.  The  effect  of  detection  * h- 
olds  has  not  been  considered  in  this  study. 

Data  from  the  USCGS  EDR  lists  for  the  first  6 months  of  1971  have  been  used  in  a study  of 
this  effect.  During  this  time  period,  1840  events  had  associated  network  m^  determinations, 
provided  by  12,049  mfe  values  from  a total  of  149  stations.  Many  of  these  stations  reported  so 
infrequently  that  they  could  be  rejected,  and,  in  fact,  the  52  stations  which  reported  more  than 
50  times  during  this  time  period  contributed  91.2  percent  of  all  m^  measurements. 

Figure  1-5  shows  the  distribution  of  USCGS  mfe  values  compared  with  that  for  mfe  values 
from  the  reduced  network  of  only  52  stations,  computed  as  the  average  of  all  magnitudes  from 
these  stations  for  each  event.  It  can  be  seen  that  the  neglect  of  the  remaining  97  stations  has 
had  little  effect  upon  either  the  shape  of  this  curve  or  the  total  number  of  events  (1825  vs  1840). 

For  each  of  these  52  stations,  biases  (station  m^  — 52-station  network  m^)  have  been  com- 
puted, and  are  listed  in  Table  1-2.  Only  events  for  which  12  or  more  stations  reported  m^  values 
have  been  used,  and  it  is  believed  that  this  limitation  ensures  that  the  average  mfe  is  not  affected 
by  a small  number  of  stations  with  large  biases.  For  each  individual  station  the  distribution  of 
the  biases  is  close  to  normal,  with  means  ranging  from  -*0.36  to  +0.31  m^  units  and  standard 
deviations  of  0.26  to  0.43.  The  stations  reporting  most  frequently  were  found  to  have  generally 
negative  biases,  which  implies  of  course  that  the  network  m^  for  the  smaller  events  for  which 
only  these  stations  report  will  be  artificially  lowered.  The  mean  biases  obtained  were  (rounded 
to  the  nearest  0.1  magnitude  unit)  applied  as  corrections  to  the  station  m^  values,  and  the  net- 
work mb  for  all  the  stations  used  was  recomputed. 

Figure  1-6  shows  the  distribution  of  m^  values  obtained  from  this  52 -station  network  with 
and  without  application  of  these  bias  corrections.  It  can  be  seen  that  the  number  of  events  at 
>4.5  has  been  reduced  as  expected,  but,  more  surprisingly,  that  the  events  have  become  re- 
distributed such  that  there  is  a substantial  increase  in  the  number  of  events  at  higher  magnitudes. 
This  difference  persists  up  to  mfe  ~ 5.7,  and  the  application  of  bias  corrections  has  increased 
the  number  of  events  at  this  magnitude  from  25  to  3 8. 

It  is  also  of  interest  to  consider  the  effect  of  removing  certain  stations  which  report  fre- 
quently. Of  the  52  stations  used,  three  (BMO,  UBO,  EUR)  in  the  western  U.S.  contributed 
23  percent  of  all  m^  observations.  Network  m^  values  have  been  recomputed,  without  applica- 
tion of  bias  corrections,  excluding  these  three  stations  only,  and  the  magnitude-frequency  curve 
for  these  is  also  shown  in  Fig.  1-6.  For  many  events,  only  these  stations  reported  m^  values 
and,  consequently,  the  total  number  of  events  has  been  reduced  to  1587,  this  difference  being 
mainly  accounted  for  by  a decrease  in  the  number  of  events  at  lower  magnitudes.  It  is,  however, 
also  evident  that  the  biases  at  these  stations  are  sufficient  to  reduce  network  mfe  values  up  to 
quite  high  magnitudes.  The  opening  or  closing  of  such  stations  can  clearly  cause  substantial 
temporal  variations  in  the  shape  of  the  magnitude-frequency  curve  up  to  quite  high  magnitude 
levels.  Furthermore,  even  if  world  network  coverage  were  ideal,  substantial  differences  in 
magnitude -frequency  curves  for  various  regions  could  be  caused  by  bias  effects. 
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TABLE  1-2 

MEAN  BIASES  OF  52  STATIONS  USED  IN  THIS  STUDY, 
WITH  RESPECT  TO  AVERAGE  OF  52-STATION  NETWORK 

Station 

Mean 

Bias 

Standard 

Deviation 

Station 

Mean 

Bias 

Standard 

Deviation 

ADE 

0.  138 

0.311 

KTG 

0.  055 

0. 247 

ALE 

0.044 

0. 239 

LAO 

0.  006 

0.369 

ALQ 

-0.011 

0.301 

LON 

-0. 291 

0. 361 

ASP 

-0. 068 

0.401 

LOR 

0.052 

0.687 

BMO 

-0.252 

0. 322 

LPB 

-0.  035 

0. 426 

BNG 

-0. 297 

0. 352 

LPS 

-0.  081 

0.338 

BNS 

0.068 

0.292 

MOX 

-0.  085 

0.293 

BUL 

-0.  067 

0.283 

NEW 

-0. 114 

0.305 

CAR 

0.  045 

0.356 

NOR 

-0.218 

0.324 

CLL 

0.  062 

0.255 

NUR 

0.  042 

0.313 

COL 

-0.044 

0. 368 

OIS 

0.  035 

0.341 

COP 

0.218 

0.280 

PMG 

0.  195 

0. 430 

CPO 

-0.  030 

0. 294 

PMR 

-0.055 

0.352 

DUG 

-0. 166 

0. 239 

PNT 

0.  123 

0. 267 

EDM 

0.308 

0. 253 

POO 

0.  124 

0.383 

EUR 

-0.358 

0.375 

PRE 

-0.  106 

0.336 

FLO 

0.288 

0. 265 

QUE 

0.  119 

0.405 

FUR 

-0.  028 

0.322 

RAB 

0.  120 

0.383 

GOL 

-0.  101 

0. 361 

RES 

0.  088 

0.316 

GRF 

0.205 

0. 334 

SJG 

0.195 

0.326 

HYB 

0.305 

0. 338 

TOL 

0.289 

0.326 

INK 

0. 139 

0. 259 

TSK 

-0. 332 

0. 407 

JCT 

0.016 

0.282 

TUC 

-0.  129 

0. 259 

KBL 

0.  067 

0. 324 

TUL 

0.207 

0.334 

KEV 

-0.024 

0.286 

UBO 

-0.  150 

0. 335 

KJN 

0.  009 

0. 293 

WIN 

-0.091 

0.317 

10 


There  are  several  possible  causes  for  the  large  but  consistent  biases  observed.  It  is  dif- 
ficult to  separate  actual  station  bias  from  the  effects  of  station  detection  thresholds,  but  it  is 
considered  that  the  values  obtained  represent  mainly  actual  bias  effects.  The  only  parts  of  the 
world  where  the  coverage  of  the  stations  used  is  sufficiently  dense  for  regional  bias  effects  to 
emerge  are  the  continental  U.S.  and  western  Europe.  Figure  1-7  shows  the  variation  of  station 
bias  across  the  U.S.,  and  it  is  clear  that  stations  in  the  basin  and  range  province  have  highly 
negative  biases  compared  with  those  further  east.  This  is  consistent  with  observations  of  the 

g 

variation  of  Q within  the  U.S.  All  stations  in  western  Europe  had  biases  of  0.0  to  0.3  m^  units, 
three  stations  in  Finland  (Baltic  Shield)  having  zero  bias.  It  is  remarkable  that  such  a poor 
measure  of  amplitude  as  can  reveal  such  differences.  Of  course,  it  is  possible  that  such 
biases  may  be  highly  azimuth-dependent  and,  in  particular,  that  certain  source-receiver  paths, 
such  as  those  traveling  down  lithospheric  slabs,  will  be  highly  perturbed.  We  propose  to  exam- 
ine this  effect  using  a larger  set  of  data  so  that  more  complete  world  coverage  can  be  attained. 

R.G.  North 

E.  mb  ESTIMATION  FOR  LARGE  EVENTS  IN  THE  PDE  CATALOGUE 

The  shape  of  the  frequency-magnitude  curve  at  its  upper  end  is  of  importance  for  several 
3 9 

reasons,  ' including  the  fact  that  this  places  a constraint  on  certain  theoretical  models  of  the 
seismic  source.  In  the  course  of  another  study  (see  Sec.  D above),  it  was  noticed  that  events  of 
very  large  m^  (>6.5)  were  not  calculated  in  a manner  consistent  with  that  used  for  smaller  events. 
In  particular,  the  PDE  magnitude  of  an  event  in  the  New  Guinea  region  (7.3)  was  calculated  as 
the  average  of  the  three  highest  m^  observations,  neglecting  eleven  which  were  less  than  m^  = 

7.0.  The  mb  observations  contributing  to  these  large  m^  values  have  been  abstracted  from  the 
earthquake  data  reports  (EDR)  of  the  PDE.  The  criteria  upon  which  PDE  m^  values  are  calcu- 
lated are  given  in  the  introduction  to  each  EDR  as:  rejection  of  any  m^  observation  for  PkP 
arrivals,  P arrivals  at  distances  of  less  than  5°,  those  associated  with  P travel-time  readings 
having  residuals  greater  than  10  sec,  and  any  differing  significantly  from  the  mean  value.  The 
rejection  of  certain  m^  observations  is,  at  least  for  the  larger  events  considered  here,  incon- 
sistent with  these  restrictions. 

We  have  recomputed  the  magnitudes  of  all  events  listed  as  having  m^  >6.5  during  the  period 
1965-1970  (inclusive),  using  the  EDR  station  magnitudes  and  our  interpretation  of  the  criteria 
listed  above.  In  every  case,  this  resulted  in  a reduction  in  the  calculated  event  magnitude.  In 
particular,  all  listed  m^  values  of  7.0  and  greater  have  been  reduced  to  below  7.0.  Since  explo- 
sions contribute  significantly  to  events  at  these  high  magnitudes,  these  too  have  been  removed 
from  the  data. 

This  revision  of  the  largest  m^  values  has  a significant  effect  on  the  shape  of  the  frequency- 
magnitude  relationship.  Figure  1-8  shows  the  relationship  for  shallow  events  in  the  PDE  listing 
for  the  period  1965-1970  (inclusive).  There  is  a suggestion  that  the  curve  becomes  linear  in  the 
mb  range  greater  than  6.5.  The  shape  of  the  curve  after  revision  of  the  large  m^  values  is  shown 
in  Fig.  1-9.  For  comparison,  this  figure  also  includes  the  same  curve  constructed  from  data  in 
the  Bulletin  of  the  International  Seismological  Center  (ISC). 

The  two  data  sets  have  a remarkably  similar  shape,  except  for  m^  < 5.0  since  the  ISC  Cata- 
logue has  a somewhat  lower  detection  threshold.  In  general,  the  ISC  Catalogue  contains  fewer 
events  with  a listed  magnitude,  since  it  requires  a minimum  of  three  station  magnitudes  before 
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an  event  magnitude  is  computed.  At  the  upper  end  of  the  curves,  the  agreement  is  excellent. 

Presumably,  the  ISC  criteria  are  very  similar  to  those  adopted  here. 

The  data  in  this  study  were  limited  to  the  period  1965-1970  since  the  magnetic-tape  version 

of  the  ISC  Bulletin  is  only  available  to  us  through  the  end  of  1970.  PDE  m^  values  for  large 

events  after  1971  appear  to  be  consistent  with  the  stated  criteria.  The  extent  to  which  these 

discrepancies  exist  at  m^^  6.5  is  not  known,  although  the  general  similarity  in  shape  of  the  ISC 

and  PDE  data  (Fig.  1-9)  suggests  that  it  is  likely  to  be  small. 

We  conclude  that  events  with  m,  > 6.8  or  6.9  occur  very  rarely,  if  At  all.  This  is  consis- 

d ^ 

tent  with  suggestions  in  earlier  studies.  R G 

M.A.  Chinnery 


F.  AMPLITUDES  AND  PERIODS  OF  THE  1965  RAT  ISLAND  SEQUENCE 

Amplitudes  and  periods  of  P-wave  first  arrivals  as  reported  in  the  Earthquake  Data  Reports 
(EDR)  for  over  300  events  of  the  Rat  Island  sequences  of  February  to  April  1965  have  been  ana- 
lyzed. Four  Vela  array  stations  — CPO  (Tennessee),  TFO  (Arizona),  WMO  (Oklahoma),  and 
UBO  (Utah)  — at  which  all  these  events  were  recorded  and  assigned  P-wave  arrival  amplitudes 
and  period  have  been  used.  These  cover  source  azimuths  of  55°  to  75°  and  are  at  distances 
of  40°  to  65°.  The  variation  in  both  azimuth  and  take-off  angles  of  the  rays  as  they  leave  the 
source  is  thus  fairly  small,  and  it  is  hoped  that  source  radiation  pattern  effects  are  negligible. 

An  immediately  noticeable  feature  of  the  data  was  that  the  reported  period  T of  the  first 
arrivals  differed  appreciably  both  at  each  individual  station  and  between  stations.  Figure  1-10 
shows  the  histograms  of  these  periods  at  each  station.  The  periods  range  from  0.3  to  2.2  sec, 
but  their  distribution  shows  large  variations  from  station  to  station.  The  most  frequently  re- 
ported periods  for  each  station  are  0.6  sec  (CPO),  0.8  sec  (TFO,  WMO),  and  0.9  sec  (UBO). 
There  is  remarkably  little  overlap  between  the  distribution  of  periods  at  CPO  and  those  at  the 
other  three  stations. 

Figure  1-11  shows  the  variation  of  log(A/T),  where  A is  the  amplitude  of  the  first  arrival, 
as  a function  of  reported  period  T for  each  station.  Each  point  is  the  mean  value  of  log(A/T), 
and  the  bar  indicates  its  associated  standard  deviation.  The  parameter  log(A/T),  with  an  addi- 
tional distance- depth  correction,  is  simply  the  body  wave  magnitude  m^.  This  figure  thus  shows 
the  variation  of  dominant  period  T with  magnitude.  In  all  cases,  T increases  quite  sharply 
with  log(A/T). 

Least-squares  fits  have  been  made  to  these  data,  and  the  values  of  a and  b where 

log  (A/T)  = a + bT  (1-23) 

and  the  correlation  coefficient  R are  given  in  Table  1-3.  Although  R varies  from  0.41  to  0.56 
reflecting  the  scatter  in  the  data,  it  can  be  seen  from  Fig.  1-11  that  these  lines  fit  the  data  quite 
well. 

The  intercept  a =0.179  for  CPO  is  ~0.2  to  0.3  m^  units  higher  than  the  other  station  inter- 
cepts, and  the  slope  b = 0.736  is  slightly  lower  than  the  slopes  of  the  other  stations.  In  Fig.  1-12, 
the  mean  values  of  log(A/T)  are  plotted  vs  T for  all  four  arrays.  Also  shown  is  the  attenuation 
effect  of  two  average  Q values  over  a distance  of  50°.  The  attenuation  is  given  by 


A(f)  = e 


-7rft* 
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TABLE  1-3 

LEAST-SQUARES  FIT  OF  DATA  TO  EQ.  (1-23)  AT  EACH  STATION. 
R IS  THE  CORRELATION  COEFFICIENT. 

Station 

a 

b 

R 

CPO 

0.179 

0. 736 

0.562 

TFO 

-0. 192 

0. 857 

0. 408 

UBO 

-0.  030 

0.857 

0.443 

WMO 

-0.003 

0.816 

0.517 

where  t*  = travel  time/Q.  The  curve  for  Q = 1330  was  calculated  using  Archambeau*s  Q 

model,*^  based  on  short-period  recordings  of  NTS  explosions.  The  t*  for  this  model  is  ~0.40, 

11 

which  is  comparable  to  other  values  measured  teleseismically. 

The  data  show  that  Q values  between  ~1330  and  665  could  explain  the  log  (A/T)  dependence 
on  period,  ignoring  seismic  scaling  effects  which  must  be  present.  Such  scaling  effects  cannot 
be  determined  independently  of  Q from  these  data.  We  plan  to  use  scaled  source  functions  from 
presumed  U.S.  explosions  on  nearby  Amchitka  Island  to  determine  Q more  accurately  for  these 
particular  paths.  This  information  will  enable  us  to  estimate  seismic  scaling  effects  for  this 
Rat  Island  sequence. 

These  data  raise  a certain  dilemma  in  that  apparent  mfe  at  these  stations  is  strongly  fre- 
quency dependent.  Thus,  reliable  mb  determinations  require  adequate  knowledge  of  Q along 

the  ray  paths.  c . W.  Frasier 

R.G.  North 
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Fig.  1-3.  Plot  of  conventional  against  ma^'mum-likelihood  m*,  estimated 
using  detection  parameters  listed  in  Table  1-1.  for  a set  of  72  events  located 
south  of  Japan  and  observed  at  a network  of  15  stations. 


Fig.  1-4.  Distribution  of  log(A/T)  for  228  events  from  Japan 
observed  at  LASA.  Smooth  curve  is  maximum-likelihood  es- 
timate. center  step  curve  is  observed  distribution,  and  others 
are  upper  and  lower  95-percent  Kolmogorov-Smirnov  confi- 
dence bands. 
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Fig.  1-5.  Magnitude-frequency  curves  for  January  through  June  1971, 
for  magnitudes  from  PDE  Catalogue,  and  from  network  of  52  stations 
uncorrected  for  bias. 


Fig.  1-6.  Magnitude-frequency  curves  for  January  through  June  1971, 
magnitudes  from  52-station  network:  (a)  uncorrected;  (b)  corrected 

for  bias;  and  (c)  uncorrected,  stations  UBO,  BMO,  EUR  removed. 
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Fig.  1-7.  Variation  of  mean  station  bias  in  across  U.S.  Dashed  line 
separates  negative  from  positive  values. 
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Fig.  1-8.  Frequency-magnitude  curve 
for  1965-1970  (inclusive)  using  values 
of  mb  listed  in  PDE  Catalogue. 
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Fig.  1-9.  Frequency-magnitude  curves 
for  1965-1970  (inclusive)  for  ISC  Cata- 
logue, and  for  PDE  Catalogue  using  re- 
vised values  of  m^  for  those  events  with 
mb  > 6.5. 
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Fig.  1-10.  Histograms  of  reported  periods  of  P-wave  arrivals  at  Vela  array  stations 
CPO,  TFO,  UBO,  and  WMO  for  302  events  in  1965  Hat  Island  sequence. 
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Fig.I-11.  Variation  of  log  (A/T)  with  period  T for  Rat  Island  events  at  four  stations 
considered.  Each  point  gives  mean  value  of  log(A/T)  at  each  period,  and  bar  indicates 
its  associated  standard  deviation. 
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Fig.  1-12.  Mean  values  of  log  (A/T)  as  a function  of  period  (from  Fig.  1-11) 
at  CPO,  TFO,  UBO,  and  WMO.  Also  shown  is  variation  of  attenuation  [ex- 
pressed as  that  of  log  (A/T)]  with  period  for  mean  path  Q values  of  1330, 
and  665. 
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II.  EARTH  STRUCTURE  AND  SCATTERING 


A.  ATTENUATION  OF  SHEAR  ENERGY  IN  THE  MANTLE 
FROM  NORMAL  MODE  ANALYSIS 

Measurements  of  Q for  a large  number  of  modes  of  free  oscillation  could  provide  the  data 

necessary  to  constrain  models  of  Q vs  depth  in  the  earth.  However,  until  now,  Q values  had 

1 2 

been  obtained  for  only  a relatively  small  number  of  modes  (e.g..  Smith,  Dratler  et  al.  ).  Thus, 

3 

meaningful  results  could  not  be  expected  from  inversion  of  these  poor-quality  data. 

A recent  advance  in  the  study  of  normal  modes  is  the  technique  of  "phase  equalization," 

4 

first  used  by  Mendiguren.  The  improved  phase -equalization  technique  used  by  Gilbert  and 
Dziewonski  (called  "stacking")  resulted  in  the  extraction  of  multiplet  shapes  for  hundreds  of 

5 

modes  that  had  not  been  detected  before.  It  was  expected  then  that  "stacking"  might  also  make 
it  possible  to  measure  Q from  the  multiplet  shapes,  and  thus  greatly  increase  the  size  of  the 
data  set  of  Q values  for  normal  modes.  Numerous  examples  of  multiplet  shapes  and  a detailed 

5 

description  of  "stacking"  and  "stripping"  techniques  can  be  found  in  Gilbert  and  Dziewonski. 

A good  "stack"  shows  a peak  corresponding  to  the  eigenperiod  of  the  particular  mode,  and  the 

shape  of  the  stack  approximates  the  spectral  shape  of  this  mode.  The  shape  depends  upon  the 

attenuation  and  can  be  used  to  measure  Q (see  Ref.  6). 

We  have  measured  Q values  for  nearly  500  modes  of  free  oscillation  by  fitting  a resonance 

th 

curve  of  the  form  of  Eq.  (H-l)  to  each  stack.  The  stack  for  the  kin  normal  mode  is  a good  approx- 
imation of  the  resonance  curve,  or  spectral  peak: 

|Ck(w)|  = [a2  + (wk -w)2)"1/2  (II-l) 

where 

ak  = V2Qk  • 

Since  we  had  reasonable  starting  values  for  and  Q^,  the  least-squares  process  converged 
after  several  iterations  to  give  the  best-fit  Q for  each  stack.  Figure  II -1  shows  examples  of 
spectral  stacks  and  the  curves  that  were  fitted  to  them.  We  note  that  Q values  obtained  in  this 
way  represent  the  minimum  Q for  each  mode  (multiplet),  since  "splitting"  due  to  the  earth's 
rotation,  ellipticity,  and  lateral  heterogeneities  broadens  the  spectral  peaks. 

For  this  study,  we  used  the  data  set  of  Gilbert  and  Dziewonski^  which  consists  of  spectra 
from  211  recordings  of  the  1963  Peru-Bolivia  and  1970  Colombia  earthquakes;  individual  record- 
ings have  an  average  length  of  about  17  hr. 

5 

Although  the  center  frequencies  could  be  estimated  for  over  800  modes,  only  477  stacks 
were  of  sufficient  quality  to  attempt  the  measurement  of  Q.  From  this  data  set,  we  consider 
the  most  reliable  Q measurements  and  coherence  of  the  Q values  obtained  for  adjacent  modes 
to  be  for  197  modes  - based  upon  the  appearance  of  the  stacks.  This  includes  49  toroidal  modes 
(only  6 of  the  fundamental  mode,  but  overtones  up  to  ^T^)  and  148  spheroidal  (with  40  QS^, 

13  < i < 66).  The  overall  period  range  of  our  data  is  from  about  500  to  80  sec. 

The  measured  Q values  for  the  fundamental  spheroidal  mode,  qS^,  agree  well  with  pub- 
lished values  determined  by  other  methods;  our  measurements  fall  within  the  scatter  of  other 
data  for  spheroidal  modes  and  long-period  Rayleigh  waves  (Fig.  II-2).  We  expect  that  the  results 
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TABLE  ll-l 

OBSERVED  Q VALUES  COMPARED  WITH  THOSE  COMPUTED 
BY  TWO  DIFFERENT  Q MODELS  OF  THE  MANTLEt 

Mode 

Q Observed 

Q Model 
(this  report) 

LMS  Q-Model 
(Smith^) 

oso 

7500  (Ref.  12) 

6402 

8816 

2S0 

704  (Ref. 2) 

1243 

3156 

4S0 

1173  (Ref. 2) 
790  (Ref.  1 ) 

1032 

2196 

5S0 

938  (Ref. 2) 
1570  (Ref.  1 ) 

1011 

2047 

6S1 

613  (Ref. 2) 

619 

1482 

8S, 

420  (Ref.  2) 

920 

2028 

11S1 

1341 

1072 

2289 

1 1 S4 

652 

669 

1426 

13S1 

1573 

940 

1942 

25S5 

791 

1010 

1579 

t Note  that  model  MLMS"  predicts  values  higher  than  observed. 
Refer  to  Fig.  11-2  for  the  parameters  of  these  models.  Numbers 
in  parentheses  following  values  of  Q observed  indicate  sources 
of  data  other  than  this  study. 
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presented  here  may  better  represent  the  average  attenuation  properties  of  the  earth.  Each  of 

our  Q measurements  is  based  upon  many  seismic  recordings,  while  previously  reported  values 

are  based  on  the  information  in  single  records. 

Several  inversions  were  performed  to  derive  Q models  of  the  mantle.  The  present  set  of 

data,  augmented  by  Q values  for  26  long-period  modes  listed  by  Smith,1  was  inverted  to  obtain 

a two-layer  model  of  the  mantle.  The  result  for  these  233  modes  was  = 120  in  the  upper 

670  km  of  the  mantle,  and  = 240  in  the  lower  2200  km.  The  variance  of  l/Q  was  such  that 

within  one  standard  deviation  of  l/Q,  Q could  vary  from  107  to  133  in  the  top  layer  and  from 

200  to  300  in  the  lower  layer.  Using  our  data  alone,  or  our  data  without  a correction  for  record 

truncation,  did  not  make  a significant  difference  in  the  result  of  the  inversion.  Observations 

and  the  theoretical  results  for  our  Q model  are  shown  in  Figs.  H-2  and  n-3  for  the  modes  0S. , 

(X  u * 

l^i'  2^i ' oTf’  l^i*  Values  f°r  a number  of  high-Q  overtones  are  listed  in  Table  II— 1 • A 
three-layer  inversion  using  the  same  data  gave  very  similar  results  for  the  lower  mantle,  and 
in  the  top  layer  (to  a depth  of  300  km)  gave  a of  86  (69  < Q < 120). 

We  have  been  assuming  that  l/Qj^,  the  dissipation  due  to  pure  dilatation,  is  zero  and  that 
the  intrinsic  Q is  independent  of  frequency.  The  former  assumption  is  justified  by  the  work  of 

7 

Birch,  and  we  tested  the  latter  by  performing  a series  of  inversions  of  our  data  within  different 
frequency  bands.  Our  data,  when  combined  with  a few  observations  of  long-period  modes  listed 
by  Smith,  cover  a period  range  from  about  1200  to  80  sec.  Since  we  have  many  observations 
of  overtones,  then,  in  any  given  frequency  band  there  are  several  modes  which  are  sensitive  to 
the  different  parts  of  the  mantle,  enough  so  that  inversions  may  be  performed  using  only  modes 
with  eigenfrequencies  in  the  same  frequency  bandwidth.  We  performed  two-layer  inversions  of 
the  data  within  16  different  frequency  bands,  each  of  a width  of  0.005  rad/sec.  A plot  of  the  re- 
sults of  these  inversions  vs  frequency,  including  the  standard -deviation  error  bars,  showed  that 
any  frequency  dependence  of  the  intrinsic  Q within  the  range  of  our  data  must  be  very  small 
(Fig.  II— 4) . 

The  fairly  low  value  of  Q for  the  lower  mantle  (200  to  300)  implied  by  our  inversions  con- 

^ 1 
trasts  with  the  higher  values  given  by  other  published  Q models  (e.  g..  Smith,  Kovach  and 

g 

Anderson  ).  Our  Q model  is  consistent  with  some  more  recent  observations  of  multiple  core 

Q 

reflections  of  ScS  which  imply  an  average  Q of  about  300  for  the  whole  mantle.  We  have  not 
yet  been  able  to  quantify  the  bias  in  our  results  due  to  splitting.  Some  preliminary  calculations 
indicate  that  this  may  not  be  as  large  as  some  have  expected.  Thus,  if  further  work  substanti- 
ates our  estimates  of  Q in  the  lower  mantle,  this  will  have  some  important  implications.  First, 
there  may  be  some  frequency  dependence  of  the  intrinsic  Q which  would  explain  the  higher 
values  obtained  for  higher  frequency  body  waves;  and  second,  the  lower  mantle  may  be  at  a high 
enough  temperature,  or  low  enough  viscosity,  to  undergo  convection. 

R.  V.  Sailor t 
A.  M.  Dziewonski 

B.  MODELING  RAYLEIGH-WAVE  PHASE  VELOCITIES  AND  AMPLITUDES  IN  EURASIA 

In  the  last  SATS,10  the  results  of  a study  of  how  well  azimuthal  deviations  of  Rayleigh-wave 
arrivals  at  NORSAR  could  be  predicted  by  ray  tracing  in  a laterally  heterogeneous  earth  were 
reported.  This  work  is  being  continued,  and  considerable  refinements  have  been  made  in  the 


t Department  of  Geological  Sciences,  Harvard  University,  Cambridge,  MA  02138. 
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grid  of  phase  velocities  across  Eurasia  used  in  the  ray  tracing.  In  particular,  different  struc- 
tural types  such  as  shield,  platform,  oceanic,  and  foldbelt  have  been  incorporated  as  well  as 

the  variation  of  crustal  thickness.  The  various  models  used  have  been  taken  from  a summary 

1 3 

of  inversions  of  observed  dispersion.  Rift-type  structures  have  not  been  included,  although 
they  are  obviously  present  in  the  Baikal  rift  zone  and  Rhinegraben,  for  example.  Much  of  the 
western  U.S.  appears  to  be  of  this  type  and  there  must  exist  large  rift  regions  in  Eurasia  which 
cannot  be  identified  from  surface  geology.  The. grid  is  specified  at  2 ° intervals  in  latitude  and 
longitude,  and  cubic  splines  are  used  during  the  ray  tracing  to  interpolate  phase  velocity  and  its 
spatial  derivatives  at  intermediate  points  and  to  ensure  their  continuity. 

We  report  the  results  of  the  ray  tracing  at  40 -sec  periods,  for  which  typical  grid  velocities 
are  3.88  (oceanic),  3.95  to  4.05  (shield),  3.35  to  3.85  (foldbelt),  and  3.5  to  4.0  (platform)  km/sec. 
Variations  in  velocity  for  a given  structure  are  caused  by  differences  in  crustal  thickness.  The 
lower  velocities  for  foldbelt  and  platform  correspond  to  the  Himalayas  and  the  Tibetan  platform, 
respectively. 

In  the  previous  SATS,1^  the  results  of  a study  of  surface  waves  from  earthquakes  in  the  Tien 

Shan  and  Pamir  mountains  were  reported.  A continuation  of  this  study  has  yielded  mechanisms 

for  these  earthquakes  from  the  phase  and  amplitude  spectra  observed  at  a network  of  stations 

1 4 

and  also  the  Rayleigh-wave  phase  velocities  for  the  event- station  paths.  These  results  place 
constraints  upon  the  velocity  grid  used,  and  it  is  also  possible  that  the  amplitude  variations  can 
be  partially  explained  by  velocity  changes  causing  focusing  or  defocusing  of  the  rays. 

Ray  tracing  has  been  carried  out  from  the  point  (39 °N,  75 °E)  in  the  Pamir  mountains.  The 
results  are  shown  in  Fig.  II-5,  a polar  plot  centered  at  the  source  and  extending  to  50°  distance. 
Also  shown  are  the  locations  of  the  stations  at  which  phase  velocity  and  amplitudes  have  been 
measured.  On  this  projection,  great-circle  paths  from  the  origin  are  straight  lines,  and  it  is 
clear  that  considerable  deviations  are  caused  by  lateral  variations  in  velocity.  Figure  II -6  shows 
on  a larger  scale  the  rays  as  they  propagate  out  to  a distance  of  15°.  The  contours  are  the  wave- 
fronts  at  phase  delays  of  100,  200,  300,  and  400  sec  due  to  propagation  effects.  In  Figs.  II-5 
and  II-6,  the  rays  leave  the  source  at  5°  increments  in  azimuth  from  0 to  355°.  The  low  phase 
velocities  in  the  Himalayas  (lower  right  of  Fig.  H-6)  can  be  seen  to  cause  considerable  delays 
and  focusing  of  the  rays.  The  extremely  large  deviations  to  the  east  in  Fig.  II-5  are  due  to  the 
fairly  rapid  transition  from  the  central  Asian  mountain  ranges  (velocity  ~ 3.7  km/sec)  to  the 
platform  region  of  North  China  (~4.0  km/sec).  The  rays  here  are  incident  upon  the  boundary 
at  low  incidence  angles  and  thus  reflection  occurs  from  the  higher  velocity  region,  causing  high 
azimuthal  deviations.  Rays  to  Northern  Scandinavia  (NW  azimuth  in  Fig.  H-5)  are  slightly  fo- 
cused by  both  near-source  velocity  variations  and  by  the  Ural  mountains;  the  latter  have  been 
assigned  a velocity  of  ~3.85  km/sec  as  opposed  to  the  platform  region  on  each  side  (4.0  km/sec). 
At  sufficiently  low  angles  of  incidence,  deviations  of  several  degrees  may  be  caused  by  rela- 
tively small  velocity  variations. 

The  ray  tracing  results  have  been  used  to  predict  phase  velocities  (distance  traveled  along 
ray/ phase  delay)  and  amplitude  variations  caused  by  nongeometric  spreading  effects  at  1 7 WWSSN 
stations  for  comparison  with  observation.  The  fit  of  the  predicted  phase  velocities  to  those  ob- 
served is  remarkably  good  [Fig.  II-7(a) ],  and  the  substantial  decrease  in  apparent  phase  velocity 
as  the  rays  sample  more  of  the  Himalayan  and  Chinese  mountain  belt  (70°  to  120°  azimuth)  is 
well  supported  by  the  data.  A slight  reduction  in  grid  velocities  in  China  should  give  even  better 
agreement.  The  variation  of  velocity  across  Europe  and  Western  Asia  (260  ° to  330°  azimuth)  is 
also  well  predicted  by  the  model. 
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The  amplitude  variation  caused  by  nongeometric  spreading  has  been  calculated  simply  by 
taking  the  ratio  of  azimuthal  spread  as  rays  leave  the  source  region  to  that  in  the  vicinity  of  the 
station.  The  observed  40 -sec  amplitudes  at  each  station  have  been  corrected  for  source  radia- 
tion pattern,  geometric  spreading,  and  an  average  Q of  200.  Predicted  and  observed  ampli- 
tudes at  each  station  are  shown  in  Fig.  II-7(b). 

The  predicted  effects  can  be  seen  to  be  quite  dramatic,  involving  a factor  of  six  in  ampli- 
tude. The  observed  variation  is  even  larger,  due  to  a combination  of  noageometric  spreading 
and  differences  in  Q.  The  large  amplitude  variation  at  European  stations  (260°  to  330°  azimuth) 
seems  likely  to  be  caused  by  lateral  velocity  variations.  The  fit  for  stations  in  Asia  is  much 
poorer,  but  it  may  well  be  that  the  revision  of  grid  velocities  in  China  required  to  fit  the  phase 
velocities  better  will  improve  the  situation. 

The  degree  to  which  the  observed  phase  velocities  can  be  predicted  is  very  encouraging  and 

it  is  hoped  that,  by  suitable  modification  to  the  grid,  the  amplitude  variations  can  be  better 

modeled.  n « 

R.  G.  North 

H.  J.  Patton 


C.  WAVENUMBER  ANALYSIS  OF  SEISMIC  CODAS 

FROM  NOVAYA  ZEMLYA  EXPLOSIONS 

In  the  previous  SATS,^  polarization  studies  of  three -component  data  from  six  Novaya 
Zemlya  explosions  showed  that  each  coda  contains  a large  number  of  impulsive  body  phases 
scattered  by  inhomogeneities  in  and  out  of  the  diametral  plane. 

We  are  now  investigating  the  wavenumber  structure  of  the  seismic  coda  using  subarray  data 
for  one  of  the  six  explosions  which  occurred  on  14  October  1969  at  7:00:06  GMT  with  an  m^  of 
6.1.  The  extended  E3  subarray  at  LASA  was  selected  to  study  the  coda  because  its  aperture  of 
20  km,  shown  in  Fig.  II-8,  gives  significantly  better  wavenumber  resolution  than  the  other  sub- 
arrays which  have  apertures  of  7 km.  The  theoretical  beam  pattern  of  the  E3  subarray  is  dis- 
played in  Fig.  II-9. 

Figures  11-10,  11-11,  and  11-12  show  3 min.  of  data  recorded  on  the  E3  subarray,  filtered  in 
three  narrow  bands  centered  at  0.7,  1.0,  and  1.3  Hz,  respectively.  At  the  bottom  of  each  figure 
is  displayed  the  unfiltered  center  trace  from  sensor  10. 

The  only  phases  after  P predicted  by  travel-time  tables  are  PcP  and  PP.  The  expected 
arrival  times  of  these  time  phases  are  shown  in  Figs.  11-10,  11-11,  and  11-12  by  dashed  lines. 
However,  it  is  clear  from  these  filtered  traces  that  other  coherent  phases  exist  which  are  at 
least  as  strong  as  the  predicted  phases. 

In  order  to  determine  the  phase  velocity  and  azimuth  of  the  coherent  arrivals  within  the  coda, 

1 5 

high -re solution  wavenumber  analysis  (described  by  Capon  ) was  applied  to  each  set  of  filtered 
traces  using  the  center  frequency  of  the  passband.  A cosine  tapered  window  10  sec  long  was  used 
to  compute  each  wavenumber  spectrum.  Successive  windows  were  made  to  overlap  by  5 sec 
so  that  no  strong  arrivals  would  be  missed  in  the  analysis. 

Figure  II-13(a-d)  shows  the  results  from  nearly  3 min.  of  coda.  Each  plotted  symbol  was 
obtained  from  high- resolution  analysis  of  a 10-sec  window  of  data.  The  time  at  which  each  svm- 
bol  is  plotted  indicates  the  center  of  the  window. 

Figures  II-13(a)  and  (b)  display  the  phase  velocity  and  azimuth,  respectively,  of  the  wave- 
number  peak,  if  at  least  2 dB  higher  than  the  background  level  on  the  plot.  Figure  II-l  3(c)  shows 
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the  amount  of  power  that  the  wavenumber  peak  exceeds  the  background  level,  and  in  Fig.  II- 13(d) 
the  average  power  of  all  25  sensors  is  plotted  for  each  window. 

The  phases  P,  PcP,  and  PP  produce  peaks  in  the  wavenumber  plots,  but  there  are  some 
interesting  variations  with  frequency.  The  P phase  arrives  with  a phase  velocity  of  18.5  km/ 
sec  for  all  three  frequencies;  however,  there  is  a monotonic  decrease  in  arriving  azimuth  with 
frequency  from  12°  to  0°.  The  wavenumber  plot  at  1 Hz  (Fig.  n-14)  shows  a very  sharp  peak  at 
18.7  km/sec,  which  is  a higher  velocity  than  the  expected  velocity  of  16  km/sec  for  P waves  at 
A = 60°.  This  velocity  shift  is  not  due  to  frequency  window  leakage  described  by  Srnart*^  be- 
cause it  is  reproducible  at  all  three  frequencies. 

The  PcP  phase  produces  strong  wavenumber  peaks  at  all  three  frequencies;  however,  the 
phase  velocities  do  not  agree.  At  1.3  Hz  the  peak  occurs  at  24.8  km/sec  (Fig.  11-15),  a little 
less  than  the  expected  27.8  km/sec,  whereas  at  0.7  and  1.0  Hz  the  wavenumber  peaks  are  much 
sharper  with  velocities  ~18. 5,  similar  to  the  P velocity.  Thus,  the  coda  at  this  time  is  a mix- 
ture of  P-type  arrivals  overlaying  the  PcP  phase,  making  it  difficult  for  the  wavenumber  anal- 
ysis to  extract  the  PcP  wave  unambiguously. 

The  PP  arrival  produces  strong  wavenumber  peaks  at  0.7  and  1.0  Hz  in  Fig.  II -1  3(c) , but  not 
at  1.3  Hz.  The  corresponding  peak  velocities  are  ~-10  to  13  km/sec,  agreeing  with  a predicted 
velocity  of  12.5  km/sec.  In  Fig.  II -16,  the  wavenumber  picture  for  the  PP  phase  is  shown  for 
0.7  Hz. 

In  Fig.  II-13(c),  we  see  that  other  strong  wavenumber  peaks  occur  at  various  times  in  the 
coda,  especially  for  0.7  Hz.  These  peaks  correspond  to  peaks  in  the  average  power  across  the 
array  shown  in  Fig.  II-  13(d).  Although  it  is  difficult  to  interpret  these  peaks  unambiguously,  it 
appears  that  most  of  these  unknown  phases  are  P-type  which  are  generated  near  Novaya  Zemlya 
within  a narrow  azimuthal  range  5°  to  10°  from  LASA.  Some  arrivals  have  velocities  which 
vary  from  PP-  to  PcP-type  arrivals,  due  to  scattering  far  from  LASA,  but  there  is  no  evidence 
of  backscattered  body  or  surface  waves  generated  locally  near  LASA. 

C.W.  Frasier 
M.  Yang 


D.  ADAPTIVE  DECONVOLUTION  OF  NOVAYA  ZEMLYA 
SHORT-PERIOD  DATA 

Deconvolution  or,  equivalently,  prediction  error  filtering  offers  a means  of  searching  for 
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reverberations  in  seismic  data  as  well  as  enhancing  first  motion.  Adaptive  deconvolution 
accomplishes  the  same  effects  while  making  as  few  assumptions  as  possible  about  the  statistics 
of  the  input  time  series.  The  filter  adapts  itself,  learning  as  it  proceeds  through  the  data,  and 
is  thus  able  to  follow  the  behavior  of  nonstationary  processes.  This  characteristic  is  attractive 
in  investigations  of  the  seismic  coda  because  both  near-source  and  near-receiver  reverbera- 
tions, as  well  as  scattering  from  distant  parts  of  the  travel  path,  may  be  involved. 

The  data  for  this  study  consist  of  three -component,  short-period  recordings  of  five  pre- 
sumed explosions  from  Novaya  Zemlya.  These  data  were  recorded  at  the  LASA  D2  subarray 
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and  were  used  previously  in  an  evaluation  of  polarization  filters.  The  locations  of  the  events 
are  shown  in  Fig.  H-17. 
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The  method  is  an  application  of  the  Widrow  algorithm  which  was  developed  by  Griffiths 
20 

et  al_.  for  use  in  petroleum  exploration.  It  can  be  expressed 

F(t  ± 1)  = F(t)  + H[x(t  + y)  - x(t  + y)]  X(t)  (H-2) 
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where  F is  the  vector  of  filter  coefficients,  x is  the  input,  x is  the  predicted  value  y time 
units  ahead,  and  X is  a vector  of  past  inputs.  The  constant  p is  a learning- rate  constant  re- 
lated to  a dimensionless  parameter  a by: 


M-  = 


0 < a < 2 


(n-3) 


2 

where  L is  the  number  of  filter  points,  and  a is  the  variance  of  the  input  data  over  the  last 
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filter  length.  Griffiths  et  aL  give  a detailed  discussion  of  how  to  choose  a and  its  effect  on 

the  filter  response.  The  value  a = 0. 5 was  used  in  this  study. 

The  only  other  two  parameters  in  the  algorithm  are  the  filter  length  L and  the  prediction 

span  y.  These  were  determined  experimentally.  Using  the  vertical  component  of  Event  1 and 

varying  the  filter  length  showed  little  change  in  the  corresponding  outputs  beyond  L = 20  points 

or  2 sec.  Furthermore,  the  prediction  error  grew  rapidly  for  any  prediction  span  beyond  1 unit 

or  0.1  sec.  Consequently,  the  values  L = 20  and  y = 1 were  used  to  produce  the  results  given 

below.  The  usual  interpretation  of  these  choices2^  is  that  the  basic  wavelet  being  deconvolved 

is  less  than  1 sec  long  and  that  the  reverberation  time  is  on  the  order  of  0.5  sec. 
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Following  Griffiths  et  a L,  the  data  were  filtered  with  the  filter  adapting  first  in  a backward 

direction  and  then,  using  the  final  filter  as  a starting  value,  in  a forward  direction.  Both  results 

were  stacked  to  produce  the  output.  The  reason  for  this  approach  is  that,  by  the  time  the  filter 

is  ready  to  go  forward  through  the  data,  it  has  already  learned  much  from  the  backward  pass. 

Some  of  the  results  for  Events  1 and  2 are  shown  in  Figs.  11-18  and  n-19.  The  principal 

feature  in  Fig.  n-18,  which  shows  the  vertical  component  plotted  in  an  expanded  time  scale,  is 

a large  negative  spike  occurring  1.8  sec  after  the  P arrival.  This  phase  was  seen  in  four  of 

the  five  events  processed,  but  only  on  the  vertical  components.  It  corresponds  to  a LASA 
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P wave  reverberation  arriving  at  approximately  1.7  sec  predicted  by  Frasier  from  the  USGS3 
crustal  model.  Unfortunately,  events  with  larger  magnitudes  recorded  at  LASA  exhibit  severe 
clipping,  which  renders  the  polarity  of  this  phase  unobservable.  Figure  n-19  shows  all  three 
components  and  their  filtered  counterparts  from  Event  2.  As  can  be  seen,  adaptive  deconvolu- 
tion removes  considerable  coda  energy  arriving  later  than  10  sec  after  the  P wave.  The  nega- 
tive spike  at  1.8  sec  can  also  be  seen  in  the  filtered  vertical  component  of  this  event. 

The  conclusions  of  this  study  are  that  adaptive  deconvolution  is  selective  enough  to  enhance 
receiver  reverberations  and,  at  the  same  time,  remove  considerable  energy  from  the  short - 
period  coda.  An  obvious  implication  of  the  latter  conclusion  is  that  much  of  this  coda  energy  is 
predictable  in  terms  of  that  which  arrives  one  filter  length  (2  sec)  before  it,  and  that  this  situa- 
tion does  not  change  faster  than  the  filter  can  adapt.  Thus,  much  of  the  short -period  coda  energy 

in  these  seismograms  appears  to  be  correlated  enough  to  be  deconvolved  out,  which  tends  to 
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support  the  earlier  conclusion  by  Greenfield  that  it  is  generated  by  scattering  in  the  source 
reglon>  D.W.  McCowan 


E.  REFLECTION  OF  PKPPKP  FROM  THE  600-km  DISCONTINUITY 

With  the  use  of  array  beamforming  and  velocity  spectral  analysis  (VESPA),  we  are  under- 
taking an  investigation  and  study  of  reflections  from  the  600-km  discontinuity.  Given  a sufficient 
number  of  these  reflections  from  a large  population  of  events  within  a given  area,  it  should  be 
possible  to  construct  a model  of  the  nature  of  the  reflecting  horizon.  The  point  of  reflection  will 
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be  located  about  240°  from  the  epicenter  of  the  event  along  the  array-event  azimuth.  For  LASA, 
for  example,  the  best  events  to  use  in  this  experiment  would  lie  about  76°  from  LASA  along  any 
azimuth.  The  resulting  PKPPKP  reflection  would  therefore  occur  about  218°  from  the  event 
along  the  same  azimuth.  Because  of  the  location  of  LASA  with  respect  to  the  seismic  zones, 
most  of  these  reflections  will  occur  in  the  ocean  basins  off  the  coast  of  Antarctica.  By  selecting 
events  in  the  Japan-Kuril  region,  it  should  be  possible  to  map  the  reflecting  600-km  boundary 
across  the  plate  boundary  of  the  Atlantic-Indian-Pacific  Basin  regions. 

Figure  11-20  shows  an  example  of  a 600-km  reflection  arriving  about  2.5  min.  before 
PKPPKP  from  a large  Novaya  Zemlya  explosion.  The  apparent  velocity  of  the  600-km  dis- 
continuity reflection  is  about  2.1  sec/deg,  while  the  apparent  velocity  of  the  PKPPKP  surface 
reflection  is  at  about  2.4  sec/deg.  The  reflection  point  for  this  event  lies  in  the  Southeast  Pa- 
cific Basin,  off  the  coast  of  Antarctica.  Because  of  the  relative  success  of  the  velocity  spectral 
analysis  process  in  detecting  the  signals  from  the  reflecting  boundary,  we  plan  to  process  a 
large  number  of  events  from  the  Kuril-Japan  region  to  investigate  the  nature  of  the  reflective 


boundary  within  a small  area. 


R.  M.  Sheppard 


F.  LATERAL  STRUCTURE  FROM  SPATIAL  CORRELATION  FUNCTIONS 
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As  part  of  the  U.S. A.- U.S.S.R.  joint  program  on  earthquake  prediction,  previous  work  on 
lateral  structure  has  been  extended  to  include  data  from  seismic  networks  in  California  and 
Garm,  U.S.S.R.  In  order  to  standardize  analysis  of  scattered  signals,  a set  of  earthquakes  in 
the  Japan  and  Kamchatka  regions  with  m^  >5.2  are  being  studied  by  determining  the  spatial 
correlation  functions  of  the  P waveforms,  of  the  time  residuals,  and  of  the  log-amplitude  re- 
siduals. These  correlation  functions  are  calculated  as  functions  of  station  separation  and  of 
distance  between  events.  Further,  correlation  functions  have  been  defined  for  distances  mea- 
sured parallel  to  the  ray  path  and  for  distances  perpendicular  to  the  ray  path;  these  functions 
have  been  termed,  respectively,  the  parallel  and  transverse  correlation  functions.  A typical 
set  of  correlation  functions  for  LASA  is  shown  in  Figs,  n-21  and  n-22.  Theoretical  considera- 
tions predict  that,  if  scattering  is  due  to  isotropic  random  inhomogeneities,  the  parallel  cor- 
relation function  will  decrease  more  slowly  than  the  transverse  correlation  coefficient.  This 
is  not  evident  in  Figs,  n-21  and  n-22,  indicating  that  the  assumption  of  isotropic  random  in- 
homogeneities is  unjustified.  Similar  results  were  obtained  for  the  preliminary  data  sets  from 
the  other  networks.  Currently,  the  remaining  correlation  functions  for  all  networks  are  being 
calculated.  j.  Scheimer 


G.  SEISMIC  EXPERIMENT  IN  TRANSMISSION  HOLOGRAPHY 

Successful  ultrasonic  model  experiments  suggest  a seismic  holography  experiment  utilizing 
NORSAR  data.  For  this  experiment,  two  events  have  been  chosen  at  moderate  distances  (77.2°, 
88.0°)  which  had  particularly  sharp  and  clear  P arrivals.  Both  events  have  nearly  the  same 
azimuth  as  seen  from  NORSAR,  but  one  event  (in  the  Japan  region,  30  October  1971)  is  rela- 
tively deep  (393  km)  while  the  other  (in  the  Kamchatka  region,  11  August  1972)  is  relatively 
shallow  (33  km). 

The  preliminary  processing  of  these  events  consists  of  filtering  records  from  the  central 
seismometers  of  each  subarray  in  three  frequency  bands:  0.5  to  0.75  Hz  (Band  1),  0.75  to  1.1  Hz 
(Band  2),  and  1.1  to  1.6  Hz  (Band  3).  For  each  of  these  filtered  records,  amplitude  and  phase 
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measurements  were  tabulated.  Prior  to  applying  the  processing  techniques  developed  in  the 
model  experiments,  the  data  were  examined  for  information  they  might  yield  on  characteristics 
of  the  medium,  and  appropriateness  of  the  frequency  bands  chosen. 

Visualizing  each  of  the  frequency  bands  contributing  a different  " color”  to  a hologram  draws 
attention  to  the  problem  of  contrast  between  each  "color."  If  the  frequency  bands  are  incorrectly 
chosen,  the  resulting  "picture"  will  be  either  (washed  out)  low  contrast  (due  to  the  similarity  of 
adjacent  frequencies)  or  will  have  such  dissimilarity  between  adjacent  frequencies  that  the  re- 
sulting monochromatic  ’^pictures"  will  be  totally  unrelated  to  one  another,  producing  a compos- 
ite that  appears  random.  The  appropriateness  of  the  choice  of  frequency  bands  can  be  examined 
by  calculating  the  correlation  coefficient  between  the  measurements  of  variations  in  log-amplitude 
for  each  adjacent  band.  A correlation  coefficient  greater  than  about  0.95  indicates  that  the  two 
sets  of  data  are  so  alike  that  little  or  no  contrast  can  be  expected.  Alternatively,  a correlation 
coefficient  less  than  about  0.5  indicates  that  the  two  sets  of  data  are  so  dissimilar  that  they  would 
not  contribute  significantly  to  each  other  when  combined.  An  optimum  range  of  correlation  coef- 
ficients might  be  from  0.75  to  0.9.  Calculated  correlation  coefficients  [p^  = p(ln  A^: In  A^)  ] be- 
tween the  filtered  records  of  the  Kamchatka  event  were: 

P12  = 0.837  , P23  = 0.841  , p13  = 0.653  . 


Rms  variations  for  the  three  bands  were  0.371,  0.495,  and  0.640,  respectively.  Mean  square 

values  were  0.137,  0.245,  and  0.410.  These  values  for  the  correlation  coefficients  indicate  that 

the  chosen  frequency  bands  should  provide  sufficient,  but  not  excessive,  contrast. 

The  reconstruction  techniques  to  be  applied  to  these  data  have  been  tested  in  concert  with 

researchers  at  the  Institute  of  Physics  of  the  Earth  in  Moscow  by  application  to  data  from  an 

ultrasonic  model  experiment.  The  model  (Fig.  11-23)  consisted  of  an  epoxy  block  in  which  was 
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embedded  an  aluminum  cross.  The  image  was  reconstructed  by  the  use  of  an  amplitude -only 

holographic  technique  using  computer  reconstruction  rather  than  the  laser  reconstruction  tech- 

24 

nique  which  has  been  widely  used  in  the  field  of  acoustic  holography.  Figure  11-24  shows  the 
results  of  this  reconstruction  of  the  wavefront  immediately  above  the  ^-rj  plane  which  contains 
the  cross.  Since  the  source  is  continuous,  the  signal  (and  reconstructed  image)  is  heavily  con- 
taminated by  surface  waves  and  multiple  reflections.  Currently  planned  are  experiments  using 
an  impulsive  source  in  order  to  fully  test  the,  numerical  techniques  prior  to  applying  them  to  the 


NORSAR  data. 


J.  Scheimer 
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Fig.  II-l.  Examples  of  spectral  stacks  for  spheroidal  modes.  Resonance  curves 
fitted  to  stacks  and  measured  Q values  are  shown.  Each  line  in  stack  is  separated 
by  7.6  x 10-6  Hz. 
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Fig.  II-2.  Theoretical  and  observed  Q values  for  fundamental  spheroidal  mode, 
0S|  . Dashed  curve  is  for  Q values  predicted  by  the  Q model  "LMS"  of  Smith* 
(Q^  = 750  from  core-mantle  boundary  to  depth  of  1000  km,  375  to  depth  of  300  km, 
and  75  in  upper  300  km).  Solid  curve  is  based  on  our  two-layer  Q model  (Q^  = 
240  in  lower  mantle,  120  in  upper  670  km),  which  was  obtained  by  inverting  Q 
measurements  for  223  modes.  Other  symbols  show  observed  Q values  from  this 
report,  Kanamori,11  and  Smith.1 
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ANGULAR  ORDER  NUMBER 


Fig.  n-3.  Theoretical  and  observed  Q vs  angular  order  number 
for  QTi#  and  1TJ?.  x's  are  observations  from  this  study, 

dashed  curved  are  values  predicted  by  model  "LMS,"  and  solid 
curves  are  values  predicted  by  our  Q model  (see  Fig.  n-2  for 
parameters  of  these  models). 


35 


1000  500  300 


PERIOD  (sec) 
3 150 


1H-M?T0I| 


eo  g 
•o  o 

40  i 

20  O 
2 


Fig.  II-4.  Results  of  inversions  using  modes  with  eigenfrequencies 
falling  within  15  separate  frequency  bands.  Error  bars  are  for  one 
standard  deviation  of  l/Q.  Dashed  lines  are  for  Q model  obtained 
by  inversion  of  223  modes.  Dots  along  bottom  of  figure  give  number 
of  modes  used  in  each  inversion. 
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Fig.  II-5.  Results  of  ray  tracing  in  grid  of  40 -sec  period  Rayleigh- wave 
phase  velocities.  Plot  is  polar  projection  centered  at  (39  °N,  74°E)  in 
Pamir  mountains,  and  extends  to  50°  distance  in  compass  directions. 
Rays  leave  source  at  5°  intervals  in  azimuth  from  0°  to  355°. 
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Fig.  H-6.  As  in  Fig.  II-5,  but  tc>  distance  of  only  15°.  Centers 
are  wavefronts  (equal  phase  delay)  at  delays  of  100,  200,  300, 
and  400  sec. 


Fig.  D-7.  Observed  and  predicted  40 -sec  Rayleigh- wave  phase  velocities 
and  amplitudes  at  17  WWSSN  stations.  Lines  joining  points  do  not  indicate 
actual  variation  of  parameters,  but  are  intended  only  to  show  matching 
characteristics  of  observation  and  theory. 
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Fig.  H-8.  Extended  E3  subarray  of  LASA. 


*43 


*34 


*73 


*43 


*73 


*54 


*63 


*84 

-22  km  — 


0.15 


tr 

o 

z 

cr 

u 

co 

35 

3 

Z 

Ui 

> 

< 

£ 


-0.15 


WAVENUMBER  EAST  (km  ') 


Fig.  II-9.  Theoretical  beam  response  of  E3  subarray.  Contour  interval  is  1 dB. 


38 


o 


Fig.  H-U.  E3  subarray  traces  of  Novaya  Zemlya  explosion,  14  October  1969,  bandpass  filtered  from  0.8  to  1.2  Hz. 


Fig.  11-12.  E3  subarray  traces  of  Novaya  Zemlya  explosion,  14  October  1969,  bandpass  filtered  from  1.1  to  1.5  Hz. 
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Fig.  H-13.  High -re solution  analysis  of  E3  subarray  data.  Each  symbol  is  result 
of  analyzing  10  sec  of  data.  Successive  windows  overlap  by  5 sec.  (a)  Shows  phase 
velocity  of  wavenumber  peak;  (b)  shows  azimuth  of  wavenumber  peak;  (c)  displays 
peak  power  above  background  on  25  sensors;  and  (d)  displays  average  power  of  all 
sensors.  All  calculations  done  at  frequencies  of  0.7,  1.0,  and  1.3  Hz. 
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Fig.  11-14.  Wavenumber  plot  of  P wave 
at  1.0  Hz.  Contour  interval  is  1 dB. 


Fig.  II-l  5.  Wavenumber  plot  of  PcP  wave 
at  1.3  Hz.  Contour  interval  is  1 dB. 


Fig.  II-16.  Wavenumber  plot  of  PP  wave 
at  0.7  Hz.  Contour  interval  is  1 dB. 
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4 12  SEP  1973  6.8 

5 27  SEP  1973  6.0 

6 27  OCT  1973  6.9 


Fig.  11-17.  Locations  (open  circles)  and  relocations  (solid  circles) 
of  five  Novaya  Zemlya  events.1 0 
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Fig.  11-18.  Original  vertical  component  and  adaptively  deconvolved  output  from  Event  1. 
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Fig.  11-19.  Original  data  and  adaptively  deconvolved  output  from  Event  2. 
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Fig.  n-20.  Array  beam  steered  for  PKPPKP  at  2.4  sec/deg  (top  trace) 
arriving  2.5  min.  after  array  beam  steered  for  PKP600PKP  at  2.1  sec/deg 
(bottom  trace).  (Gains  and  signals  are  adjusted  for  visual  comparison.) 
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Fig.  11-21.  Spatial  correlation  functions 
derived  for  LASA  from  set  of  Japan 
earthquakes. 


Fig.  11-22.  Spatial  correlation  functions 
derived  for  LASA  from  set  of  Kamchatka 
earthquakes.  In  both  cases,  parallel  and 
transverse  refer  to  station  separations 
projected  onto  axes,  respectively,  paral- 
lel and  perpendicular  to  ray  path. 
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Fig.  II -2 3.  Diagram  of  epoxy  model  used  in  ultrasonic  holography 
experiment.  £-tj  plane  lies  30 A.  above  source,  s,  and  20A  below 
x-y  plane.  Data  were  recorded  at  x-y  plane  every  l/3X  to  give 
matrix  of  6l  X 61  data  points  from  which  wavefront  immediately 
above  £ -tj  plane  was  reconstructed.  Source  is  12A.  diameter. 


Fig.  11-24.  Contour  map  of  amplitudes  of  reconstructed  wavefront  just 
above  £-t)  plane  normalized  to  arbitrary  maximum  of  2.0.  Shaded  area 
corresponds  to  position  of  cross-shaped  inhomogeneity.  Salient  features 
of  cross  are  clearly  depicted  by  reconstructed  wavefront. 
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III.  MISCELLANEOUS  STUDIES 


A.  STUDIES  ON  THE  MAXIMUM  ENTROPY  SPECTRAL  TECHNIQUE 

An  important  characterization  of  many  physical  processes  is  the  power  spectrum.  Most 

standard  techniques  used  to  estimate  this  function  weight  the  data  in  such  a way  that  the  data  or 

their  autocorrelation  are  zeroed  beyond  some  lag.  When  the  physics  of  the  process  dictates 

stationarity,  such  techniques  can  compromise  absolute  resolution  and  produce  spectral  overlap 

of  adjacent  components.  Maximum  Entropy  Spectral  Analysis,  on  the  other  hand,  estimates  the 

spectrum  without  assuming  zero  data  outside  the  observed  region.  This  technique,  using  Burgs' 

innovation  for  calculating  a consistent  prediction  error  wavelet  directly  from  observed  data,  has 

1 2 

proven  very  successful  in  estimating  spectra  in  several  seismologically  important  areas.  * 

There  are,  though,  several  ill-defined  problems  in  applying  the  technique  to  data  with  sig- 
nificant amounts  of  noise.  Most  important  is  use  of  the  proper  filter  length.  The  problem  is 

currently  being  studied  using  the  correspondence  between  the  Maximum  Entropy  Method  and 

3 4 

autoregressive  techniques  noted  by  Ulrych  and  Bishop.  Akaike  has  proposed  two  methods  for 
determining  the  Order  of  an  autoregressive  process.  Essentially,  the  error  in  predicting  the 
data  as  a function  of  filter  length  is  weighted  such  that  less -than- random  decreases  in  the  error 
produce  an  increase  in  the  weighted  error.  The  estimate  of  the  order  then  lies  at  the  minimum 
of  the  functions.  Figure  III-l  shows  a representative  error  curve  and  weighted  error  curves 
used  in  the  Maximum  Entropy  Cepstral  Technique  shown  elsewhere  in  this  report.  Here,  a 
triple  check  of  the  technique  is  available  since  the  real,  imaginary,  and  complex  log  spectra 
should  all  have  the  same  order.  Figure  III— 1 , showing  the  results  for  just  the  imaginary  log 
spectrum,  indicates  an  order  of  23.  The  real  and  complex  log  spectra  had  the  same  values. 
Computed  spectra  using  this  filter  length  gave  reasonable  results.  Unfortunately,  data  with 
larger  amounts  of  noise  do  not  give  consistent  results,  and  research  into  this  problem  is 
continuing.  T#  E.  Landers 

B.  THE  VARIATION  OF  MAXIMUM  ENTROPY  CEPSTRAL  pP  TIMES 

OVER  A NETWORK  OF  STATIONS 

In  anticipation  of  SRO  data,  the  raw  traces  from  the  short-period  center  instruments  of 
NORSAR  for  an  NTS  explosion  and  a presumed  explosion  in  Eastern  Kazakh  (EK)  have  been 
analyzed  by  the  Maximum  Entropy  Cepstral1  method  to  determine  the  variation  of  pP  delay 
times  observed  over  a network  of  stations.  Since  waveforms  rather  than  absolute  arrival  times 
or  amplitudes  comprise  the  data  set,  use  of  the  individual  stations  as  a network  is  justified 
despite  all  stations  being  approximately  equidistant  from  the  events.  Short-period  waveform 
degradation  due  to  near-receiver  scattering  is  severe  across  NORSAR,  and  the  ability  of  any 
individual  waveform  technique  to  give  consistent  results  will  be  adequately  tested.  Since  each 
path  or  near-receiver  scatterer  has  approximately  the  same  effect  on  P and  pP,  the  basic 
echo  delay  time  should  be  well  determined  even  though  signal  shapes  vary  significantly.  Un- 
fortunately, dispersed  waveforms  have  poorer  signal-to-noise  (S/N)  ratios  and  so  less-well- 
determined  cepstra. 

The  source  parameters  for  the  events  used  in  this  study  were  obtained  from  the  PDE  and 
NORSAR  Bulletins,  and  are  listed  in  Table  III— 1.  Data  from  each  center  site  of  the  C-ring  sub- 
arrays and  the  OA  subarray  were  used.  Approximately  15  sec  of  each  of  the  15  signals  for  each 
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TABLE  lll-l 

SOURCE  PARAMETERS  OF  THE  EVENTS  USED  IN  THE  STUDY 

Event 

Date 

(1974) 

Time 

Latitude 

Longitude 

Depth 

mb 

A 

(deg) 

Azimuth 

(deg) 

NTS 

EK 

20  Aug 
30  Jan 

15:00:00 

04:57:02 

37°  N 
50°N 

1 16°W 
78°E 

OG 

OG 

5.8 

5.4 

73 

38 

318 

75 

event  was  analyzed.  An  example  of  the  waveforms  is  given  in  Fig.  Ill— 2.  Whereas  the  longer- 

period  NTS  event  signals  remain  relatively  constant,  the  shorter-period  EK  event  records  at 

the  same  two  sites  show  the  strong  dispersion  of  signal  energy  that  is  often  observed.  The 

result  is  that  while  the  S/N  ratio  at  OC9  is  approximately  10/1,  at  OC13  it  is  only  2/1.  Similar, 

though  not  so  severe,  variations  also  occur  for  event  NTS. 

1 

In  the  manner  described  in  the  previous  SATS,  the  Maximum  Entropy  Cepstrum  was  com- 
puted for  each  waveform.  For  each  event,  the  log-spectrum  band  and  t*  were  held  constant  in 

4 

the  analysis  of  each  waveform.  In  addition,  Akaike’s  order  criterion  was  computed  as  an  aid 
in  determining  the  filter  length.  Cepstral  peaks  were  read  to  the  nearest  digitizing  interval, 

0.05  sec.  The  results  are  given  in  Table  m-2,  and  are  shown  in  histogram  form  in  Fig.  Ill— 3 . 


TABLE  III -2 

pP  TIMES  OBSERVED  AT  THE  NORSAR 

SUBARRAYS 

Event 

A1 

Cl 

C2 

C3 

C4 

C5 

C6 

C7 

C8 

C9 

CIO 

Cll 

C12 

C13 

C14 

NTS 

0.80 

0.70 

0.80 

0.70 

0.80 

0.85 

0.80 

0.75 

0.90 

0.90 

0. 85 

0.80 

0.75 

0.75 

0.80 

EK 

0.45 

0.70 

0.65 

0.60 

0.55 

0.75 

0.60 

0.55 

0.65 

0.55 

0.50 

0.70 

0.60 

0.55 

0.65 

The  mean  and  standard  deviation  of  event  EK  is  0.60  ± 0.08  sec,  and  for  event  NTS  it  is  0.80  ± 
0.06  sec.  The  consistency  of  these  results,  despite  the  complicating  receiver  effects  and  large 
S/N  ratios,  are  very  encouraging. 

T.  E.  Landers 

C.  EARTHQUAKE  SOURCE  MECHANISM  DISCRIMINATION 

Work  continues  on  the  study  of  discrimination  between  earthquake  source  mechanism  types. 

1 

As  reported  in  the  last  SATS,  the  mid- Atlantic  region  is  studied  because  the  earthquakes  are 
predominantly  of  two  types  (strike-slip  and  normal),  both  occurring  at  similar,  shallow  depths. 
Events  of  this  source  region  are  studied  to  obtain  a better  understanding  of  the  typical  spread 
observed  in  explosion -earthquake  discriminants,  such  as  Mg  vs  m^. 

The  separation  observed  in  Mg  vs  m^  for  the  two  mechanism  types  was  attributed  to  the 
combination  of  the  effects  of  source  geometry  and  a local  zone  of  attenuation  associated  with  the 
spreading  ridge.  Dip-slip  geometry  causes  m^  to  be  higher,  while  an  attenuating  zone  under 
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the  event  reduces  However,  the  effect  of  scattering  on  must  also  be  considered.  If  a 

significant  amount  of  the  short-period  energy  radiated  from  an  event  undergoes  scattering,  the 
lobes  of  the  radiation  pattern  will  effectively  be  smoothed  out.  As  a result,  the  m^  for  dip-slip 
events  will  not  be  increased  as  greatly  as  anticipated.  This  effect  of  scattering  is  being  studied 
to  re-examine  how  much  change  in  Mg  vs  m^  can  actually  be  attributed  simply  to  the  effect  of 
the  radiation  patterns  for  particular  source  geometries. 

Although  signal  complexity  is  also  a useful  -parameter  in  the  explosion  discrimination  prob- 
lem, it  has  not  been  found  useful  in  separating  these  two  source  mechanism  types.  This  is 
surprising  in  view  of  the  observed  source  complexity  of  shallow  strike -slip  events  (noted  by 
Filson,5  for  example).  Work  continues  on  the  modification  of  the  traditional  complexity  measure 
to  obtain  a parameter  which  does  reflect  this  observable  mechanism -dependent  signal  character. 

A.  F.  Shakal 

D.  PARAMETRIC  TRAVEL-TIME  TABLES  FOR  TELESEISMIC  PHASES 

An  important  motivation  for  construction  of  a Parametric  Earth  Model  (PEM)^  has  been 
that  its  functionals  such  as  travel  times  are,  for  a particular  branch,  smooth  functions  of  the 
epicentral  distance.  The  density  and  velocities  in  model  PEM  are  described  by  piecewise  con- 
tinuous functions  of  radius;  polynomials  up  to  the  third  order.  Therefore,  the  travel  times  for 

this  model  might  also  be  represented  in  a similar  form;  such  representation  has  been  shown  to 

7 8 

be  appropriate  in  smoothing  the  observed  P and  S travel  times.  ' Model  PEM  is  a good  fit  to 

9 

an  extensive  set  of  free  oscillation  data  as  well  as  to  several  subsets  of  the  body  wave  travel 
times:  P,  S,  SKS,  SKKS,  and  PKIKP.  It  may  be  expected,  therefore,  that  the  travel  times 
computed  for  this  model  would  be  more  representative  of  the  real  earth  than  the  Jeffreys-Bullen 
tables*^  computed  nearly  four  decades  ago.  The  advantages  of  an  analytical  representation  of 
the  travel  times  are  the  following:  (1)  the  travel  times  for  a particular  branch  are  defined 
uniquely,  there  is  no  need  for  interpolation  which  is  always  subject  to  a number  of  assumptions; 
(2)  the  functions  can  be  analytically  differentiated  to  obtain  the  slowness  parameter  and  ampli- 
tudes; (3)  computer  storage  requirements  are  significantly  reduced;  and  (4)  certain  other 
parameters  such  as  the  incidence  angles  or  bottoming  depths  can  also  be  computed  using 
analytical  formulas. 

Coefficients  of  the  polynomial  representation  of  the  teleseismic  P travel-time  branch  (rays 
bottoming  below  650-km  discontinuity)  are  shown  in  Table  m-3;  all  computations  are  performed 
for  the  continental  model  PEM-C. 

The  formula  for  evaluation  of  the  travel  time  is 

N 

1(A)  = £ anxn  (III-l) 

n=0 

where  x = (A  — B)/A. 

The  constants  A and  B are  chosen  such  that  x assumes  values  between  —1  and  1,  approxi- 
mately. This  means  that  all  the  coefficients  a^  need  to  be  specified  with  the  same  precision. 

The  values  of  the  polynomial  coefficients  are  obtained  by  the  least -squares  fit  of  p(A)  com- 
puted for  an  appropriately  selected  set  of  ray  parameters;  this  assures  that  p(A)  obtained  by 


51 


TABLE  1 II  — 3 

PARAMETRIC  REPRESENTATION  OF  THE  P TRAVEL  TIMES 
FOR  RAYS  BOTTOMING  IN  THE  LOWER  MANTLE^ 

Depth 

(km) 

A . 
min 

A 

max 

°0 

ai 

°2 

°3 

a4 

°5 

°6 

°7 

0 

17.2 

94.7 

609.00 

274.56 

-60.24 

2.55 

3.75 

-5.72 

1.11 

1.64 

20 

17.1 

94.6 

605.76 

274. 36 

-60.18 

2.56 

3.67 

-5.69 

1.13 

1.62 

35 

17.0 

94.6 

603.65 

274.18 

-60.12 

2.57 

3.62 

-5.65 

1.13 

1.61 

70 

16.8 

94.5 

599.88 

273.63 

-59.92 

2.58 

3.45 

-5.54 

1.17 

1.57 

120 

16.3 

94.3 

594.48 

272.84 

-59.63 

2.61 

3.23 

-5.35 

1.22 

1.52 

170 

15.9 

94.2 

588.93 

272. 07 

-59.36 

2.62 

3.01 

-5.17 

1.27 

1.46 

220 

15.5 

94.0 

583.39 

271.29 

-59.08 

2.63 

2.79 

-5.00 

1.30 

1.40 

270 

15.0 

93.8 

578.56 

270.39 

-58.74 

2.60 

2.57 

-4. 77 

1.33 

1.33 

320 

14.5 

93.7 

573.84 

269.46 

-58.37 

2.57 

2.36 

-4.53 

1.35 

1.27 

370 

13.9 

93.5 

569.25 

268.51 

-58.00 

2.50 

2.15 

-4.27 

1.36 

1.19 

420 

13.3 

93.3 

564.75 

267.51 

-57.61 

2.42 

1.96 

-4.02 

1.36 

1.11 

470 

12.5 

93.0 

560.68 

266.42 

-57.16 

2.31 

1.81 

-3.70 

1.33 

1.01 

520 

11.7 

92.8 

556.66 

265.31 

-56.70 

2.16 

1.67 

-3.38 

1.28 

0.92 

570 

10.8 

92.6 

552.73 

264.17 

-56.24 

2.00 

1.57 

-3.07 

1.20 

0.81 

620 

9.8 

92.3 

548.87 

262. 99 

-55.76 

1.83 

1.51 

-2.77 

1.11 

0.72 

670 

8.7 

92.1 

545.09 

261.79 

-55.27 

1.64 

1.49 

-2.48 

1.01 

0.63 

720 

14.4 

91.8 

541.94 

260. 38 

-54.70 

1.12 

1.97 

-1.74 

0.58 

0.11 

770 

16.4 

91.5 

538.86 

258.92 

-54.01 

0.69 

2.09 

-1.30 

0.47 

-0.20 

t The  constants  A and  B in  x = (A  — B)/A  are  40°  and  60°,  respectively.  For  foci  shallower 
than  670  km,  the  travel  times  for  distances  near  A^.^  will  correspond  to  secondary  arrivals. 

52 


differentiation  of  Eq.  (III-l)  will  be  accurate: 


N 

p(A)  = A"1  Yj  nanxn_1  (HI-2) 

n=l 


coefficient  a^  is  then  obtained  to  give  the  best  least-squares  to  the  travel  times. 

For  the  set  of  parameters  shown  in  Table  III-3,  the  rms  difference  between  p(A)  obtained 
using  Eq.  (Ill- 2)  and  computed  exactly  is  0.002  sec/deg  and  the  maximum  difference  is 
—0.008  sec/deg  at  94.7°;  the  rms  error  for  the  travel  times  is  0.03  sec,  and  maximum 
difference  0.07  sec  at  93°.  Clearly,  the  fit  is  satisfactory. 

Travel  times  have  also  been  calculated  for  a number  of  other  phases.  For  instance,  poly- 
nomial coefficients  for  the  surface  focus  are  listed  in  Table  III -4  for  52  phases;  the  full  set  is 
much  more  extensive,  as  the  SP  and  PS  travel  times  are  different  for  a finite  depth  focus;  also, 
multiple  surface  reflections,  trivial  for  a surface  focus,  are  not  included.  Note  that  only  the 
phases  with  a ray,  or  both  rays  in  the  case  of  converted  phases,  bottoming  in  the  lower  mantle 
are  listed.  The  discontinuity  in  the  upper  mantle  of  model  PEM  gives  rise  to  multiple  branches 
for  rays  bottoming  in  this  region;  the  best  representation  for  the  upper-mantle  travel  times  is 
still  under  consideration. 

In  addition  to  the  travel  times,  ellipticity  corrections  have  also  been  computed  using  the 

111 

approach  of  Dziewonski  and  Gilbert  ' in  which  the  ellipticity  correction  is  represented  by: 

2 

6t=  £ P2>m(eo)  ' cos  * Tm^)  ’ (HI-3) 

m=0 


The  analytical  representation  of  these  corrections  in  terms  of  Eq.  (III-l)  is  possible,  al- 
though it  is  more  complex  since  the  functions  Tm(A),  unlike  the  travel  times,  are  generally  not 
monotonic.  For  some  phases  (e.g.,  SKKKS),  rm(A)  have  four  local  extrema;  accurate  repre- 
sentation with  a low  degree  polynomial  is  not  possible  in  this  case. 

This  difficulty  may  be  circumvented  in  the  following  way.  The  most  significant  part  of  the 
ellipticity  correction  arises  due  to  a change  in  the  radius  of  the  free  surface;  this  is  particularly 
true  for  the  core  phases.  Thus,  instead  of  the  form 

N 

V|A>-  2 »nm*n  (III-4) 

n=0 


we  may  seek  the  representation: 


r (A)  = C [6  n + P0  (A)] 
m m mO  2,m 


N 

+ l 

n=0 


c x 
nm 


(III- 5) 


[see  Eqs.  (23)  to  (25)  in  Ref.  11];  Eq.  (Ill- 5)  can  easily  be  generalized  for  the  multiple  surface 
reflections.  In  Table  III-5,  we  compare  results  obtained  for  the  core  P phases  using  Eqs.  (III-4) 
and  (III- 5).  For  PKIKP  and  PKP^,  nearly  the  entire  variation  of  Tm(A)  is  described  by  the 
first  term  in  Eq.  (in- 5 ).  This  is  in  good  agreement  with  the  J-B  tables,*^  where  a constant 
value  of  f(A)  = 0.09  sec  is  suggested  for  all  the  PKP  branches.  This  value  of  f(A)  is  equivalent 
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. 

PARAMETRIC  REPRESENTATION  OF 

TABLE  1 11—4 
TRAVEL  TIMES  FOR  52 

PHASES  - 

SURFACE  FOCUS 

Phase^ 

A 

B 

A . 
min 

A 

max 

ao 

°i 

a2 

°3 

a4 

°5 

°6 

°7 

PKP 

30 

150 

114.2 

180.0 

1187.75 

46.86 

-14.88 

-8.91 

0.45 

2.42 

-0.13 

-0.50 

PKPj 

5 

150 

144.7 

155.5 

1193.01 

12.78 

-1.27 

0.14 

0. 00 

-0.01 

-0.05 

0.04 

pkp2 

15 

155 

144.7 

173.1 

1220.04 

64.33 

2.32 

-1.04 

0.41 

-1.27 

1.47 

-0.50 

PcP 

45 

45 

0.0 

94.3 

598.27 

155.80 

47.40 

-21.62 

2.33 

1.73 

-0.76 

- 

PKiKP 

55 

55 

0.0 

113.6 

1028.34 

63.20 

27.24 

-4.02 

-0. 68 

0.14 

- 

- 

PKKP 

80 

280 

207.4 

360.0 

1845.61 

128.41 

-40.56 

-26.69 

2.00 

8.85 

-0.54 

-2.30 

PKKP] 

25 

260 

235.8 

288.0 

1806.34 

71.95 

-10.09 

0.  91 

0.08 

-0. 06 

-0.47 

0.29 

pkkp2 

10 

240 

235.8 

251.9 

1744.24 

43.64 

0.82 

-0.49 

0.12 

- 

- 

- 

PKKKP 

120 

420 

300.6 

540.0 

2519.26 

183.76 

-65.22 

-33.18 

6.98 

9. 02 

-2.05 

-2.06 

PKKKP] 

50 

370 

320.1 

420.5 

2419.11 

146.47 

-23.24 

2.09 

0.24 

-0.83 

-0.85 

0.86 

pkkkp2 

5 

325 

320.1 

330.6 

2268.78 

22.10 

0.06 

-0.04 

0.01 

- 

- 

- 

PKIIKP 

120 

240 

114.7 

360.0 

1342.70 

166.81 

-65.04 

-19.52 

5.19 

1.32 

-0.89 

- 

PcPPKP 

25 

160 

135.2 

180.0 

1700.93 

44. 72 

-5.38 

-8.00 

-11.28 

-6. 00 

4.03 

3.65 

PcPPKP] 

45 

220 

178.4 

267.4 

1894.87 

189.25 

12.49 

-6.91 

-0.08 

-0.78 

5.00 

-2.83 

PKPPKP670 

70 

290 

226.3 

360.0 

2206.16 

117.59 

-30.01 

-26.41 

-3.89 

9.41 

1.57 

-2.72 

PKPPKPj670 

10 

295 

285.0 

308.6 

2222.24 

26.64 

-2.56 

0.25 

-0.01 

0. 08 

-0.13 

0.05 

PKPPKPj70 

30 

310 

285.0 

341.0 

2295.99 

130. 02 

3.71 

-2.00 

-0.41 

-0.27 

2.22 

-1.26 

t The  superscripts  "670"  denote  the  underside  reflections  from  the  670-km  discontinuity, 
traveling  through  the  inner  core. 

The  core  ph 

ases  without  i 

a subscript  refer  to  rays 

TABLE  MI-4  (Continued) 

Phase 

A 

B 

A . 
min 

A 

max 

ao 

ai 

a2 

°3 

a4 

°5 

a6 

°7 

PR670 

80 

100 

17.4 

184.2 

952.34 

579.18 

-111.03 

-3.25 

10. 37 

-6.32 

-0.19 

1.25 

s 

40 

60 

15.4 

98.9 

1106.72 

514.43 

-91.68 

-0.61 

5.33 

-6.29 

0.50 

1.32 

ScS 

50 

50 

0.0 

98.7 

1132.75 

342.67 

95.26 

-50.68 

8.37 

3.84 

-2.36 

- 

SKS 

40 

140 

104.1 

180.0 

1604.93 

64.67 

-19.89 

-13.54 

0.61 

4.47 

-0.14 

-1.17 

SKS] 

35 

110 

70.6 

144.3 

1515.73 

144.25 

-43.60 

7.64 

-4.83 

1.95 

2.31 

-1.90 

SKKS 

80 

280 

197.3 

360.0 

2278.41 

120.02 

-44.10 

-20.15 

5.10 

4.86 

-1.44 

-1.05 

SKKSj 

95 

180 

81.1 

276.9 

1974.38 

434.04 

-138.22 

19.22 

-0.43 

-6.48 

4.04 

-0.65 

SKKKS 

125 

415 

290.5 

540.0 

2944.12 

188.37 

-69.35 

-31.84 

8.46 

7.51 

-2.46 

-1.52 

SKKKSj 

160 

250 

91.6 

409.4 

2431.32 

749.21 

-246.04 

29.14 

5.00 

-13.62 

6.03 

-0.35 

ScSSKS1 

15 

140 

128.9 

156.1 

2526.59 

52.36 

-12.33 

1.93 

-0. 30 

1.93 

-3.18 

1.35 

ScSSKS2 

2 

130 

128.9 

131.6 

2486. 1 1 

13. 30 

0.65 

-0.04 

0.42 

-0.65 

-0.20 

0.47 

PS 

15 

125 

110.1 

137.8 

1857.58 

132.89 

-3.33 

-0.77 

0.17 

- 

- 

- 

PPS 

25 

150 

127.3 

176.6 

2183.70 

223.41 

-4.48 

-1.39 

0.17 

- 

- 

- 

PSS 

20 

220 

203.0 

236.7 

3344.55 

177.21 

-4.73 

-1.04 

0. 36 

- 

- 

- 

PcS 

30 

30 

0.0 

60.7 

777.26 

97. 76 

35.46 

-11.64 

-1.11 

0. 73 

- 

- 

PKS 

35 

155 

109.1 

180.0 

1410.62 

43.22 

-26.00 

-8.05 

4.69 

1.80 

-1.89 

-0.88 

PKS] 

10 

140 

130.0 

147.6 

1389.16 

27.49 

-3.79 

0.46 

-0.10 

-0.22 

-0.05 

0.26 

pks2 

5 

135 

130.0 

139.5 

1378.95 

22.00 

0.16 

-0.08 

0.07 

-0. 05 

- 

- 

PScS] 

5 

115 

110.3 

121.8 

1771.62 

27.02 

-2.19 

0.24 

0.01 

0.10 

-0.18 

0.07 

PScS2 

15 

125 

110.3 

137.6 

1861.06 

124.53 

2.23 

-1.85 

0.14 

0.91 

0.72 

-1.12 
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PARAMETERI 

ZATK 

TABLE  111-5 

ON  OF  ELLIPTIC1TY  CORRECTION  COEFFICIENTS  FOR  PKP  PHASES* 

Phase 

and  A Range 

m 

C 

m 

o 

o 

-o 

Vcl 

b2/c2 

b3/c3 

b4/c4 

b5/c5 

b6/c6 

0 

-2.10 

-1.05 

0.25 

0.21 

0.00 

-0.01 

-0.01 

-1.35 

0.00 

- 

| . 

1 _ 

PKIKP 

1 

- 

1.07 

-0.58 

-0.58 

0.09 

-0.14 

0.01 

0.13 

114.2  to  180.0 

i 

-1.27 

0.07 

-0.07 

1 

| 

- 

1 

! - 

- 

2 

-0.32 

0.53 

-0.12 

-0.02 

-0.02 

-0.07 

0.01 

— w 

-1.22 

■ • 

0.00 

~ 

L 

- 

i - 

- 

0 

_ 

-2.16 

j -0.20 

0.02 

i -0.01 

0.01 

- 

- 

-1.33 

0.00 

-0.04 

0.01 

1 

- 

i 

- 

PKP] 

1 

-• 

0.99 

-0.09 

• -0.02 

0.01 

-0.01 

I 

- 

- 

144.7  to  155.5 

-1.32 

0.00 

0.02 

-0.01 

- 

1 

1 

- 

2 ! 

- 

-0.29 

0.0? 

- 

- 

-1.33 

1 

0.00 

_ 

0.00 

-o.oi  ! 

_ 

_ 

C 

-2.06 

-0.19 

0.05 

0.00 

0.02 

0.03 

-0.03 

-1.14 

-0.02 

C.  10 

1 

-0.03 

0.01 

-0.01 

0.01 

- 

PKP2 

1 

- 

0.63 

-0.39 

-G.03 

0.01 

-0.01 

-0.02 

0.02 

144.7  to  173.1 

-1.10 

-0.02 

-0.04 

0.03 

-0.01  1 

0.03 

1 

O 

8 

- 

1 

2 

- 

-0.35 

0.05 

-0.03 

0.00 

-0.01 

-0.01 

0.01 

1 

I 

1 

-1.37 

-0.20 

-0.15 

0.02 

0.00 

0.01 

mm 



L 

t Obtained  using  Eg.  (1 11-1 ) - constants  b^,  and  Eq.  (111-5)  - 

table:  A = (A  - A . J/2-  B = (A  ♦ A . )/2. 

max  mm  v • ' mo*  min'  ’ 

constants  Cm  and  c^.  Por  this 

S7 


to  C = -1.28  in  Table  Ql-S,  with  all  other  coefficients  being  zero.  The  maximum  errors 
m 

resulting  from  this  approximation  would  be  of  the  order  of  0.1  sec  for  PK1KP  and  PKP^,  but 
for  PKP,  errors  could  amount  to  as  much  as  0.5  sec. 

A.  M.  Dziewonski 


E.  TOWARD  A MOMENT  TENSOR -ENERGY  RELATION 

The  advantages  to  expressing  seismic  radiation  fields  in  terms  of  the  moment  tensor  of  the 
applied  source12  are  well  known.  The  six  independent  quantities  in  this  tensor  provide  informa- 
tion about  the  seismic  source  in  the  most  general  form;  for  example,  it  can  describe  a source 
representing  a linear  combination  of  three  point  source  mechanisms;  explosive,  compensated 
linear  vector  dipole,  and  double  couple.13  Furthermore,  the  time  history  of  the  moment  tensor 
determines  the  phenomenology  of  the  event.14  In  this  report,  we  show  howto  solve  the  forward 
problem  of  directly  calculating  the  radiated  seismic  energy  if  the  moment  tensor  (i.e..  source 
mechanism)  of  an  event  and  the  earth  structure  are  known.  In  our  approach,  we  utilize  certain 
fundamental  properties  of  normal  modes. 

Following  Gilbert.12  we  can  express  the  vector  displacement  of  the  earth’s  surface  as  a 
sum  over  its  normal  modes  of  oscillation: 

1 -cosu)  t 

s I lM:Sn(*Vl  — “ l * (ID’ 

n FnWn 

In  this  expression,  s^  is  the  displacement  vector  for  the  n1*1  mode,  M is  the  moment  tensor  of 
the  source.  Sn(Fg)  is"the  strain  tensor  of  the  n01  mode  evaluated  at  the  source  position.  «n  is 
the  angular  frequency  of  the  nth  mode,  and  FR  is  the  normalization  factor: 

F = ^ P(F)|  s*(F)|2  dv  . 
n 

The  kinetic  energy  density  in  each  mode  will  then  be  given  by: 

I sinw  ti2 

f-n(r)  = ^ p(F)|rn(F)Wn  -y-z- 1 


(IU-7) 


(111-8) 


where  Wn  represents  the  0th  weight:  S:M.  Integrating  this  over  the  volume  of  the  earth  and 
evaluating  it  at  its  maximum  (u>nt  = »/2)  gives  the  kinetic  energy  in  each  mode: 


E 


n 


_1 

2 


F u:  1 
n n 


(HI-9) 


where  we  have  used  the  normalization  relation  (II1-7)  in  the  numerator.  Since  the  normal  modes 
are  independent,  the  total  energy  is  the  sum  of  the  maximum  kinetic  energies: 


E 


n 


M_S^/ 

Fn"n 


(111-10) 


This  relationship  can  be  generalized  for  an  event  with  an  arbitrary  time  history  if  M in 
Eq.  (iu-tO)  were  to  represent  the  spectrum  of  the  moment  rate  tensor  M(w)<  For  the  case  of 
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a ramp  increase  in  applied  moment  of  length  T,  and  the  total  moment  Mq: 


i y 

1 Mo:S„irs>l 

/sin  -§- 

2 u 

F w 2 

\ "nT 

n 

n n 

\ 2 

(III-ll) 


As  can  be  seen,  the  effect  of  changing  the  source  time  function  from  a step  function  to  a 
ramp  function  is  to  reduce  the  total  radiated  energy.  Thus,  events  with  slow  increases  in 
applied  moment  radiate  decreasing  amounts  of  high-frequency  energy.  Energy  associated  with 
the  low-frequency  modes  is  relatively  unaffected. 

To  apply  either  of  these  expressions  to  actual  events,  it  is  necessary  to  have  a complete 
catalog  of  normal  mode  frequencies  and  their  associated  strains.  A suitable  model  which  satis- 

Q 

fied  a wide  spectrum  of  standard  earth  data  is  Gilbert's  and  Dziewonski's  1066 A model  for 

density,  and  compressional  and  shear  velocities.  For  this  model,  we  have  spheroidal  modes 

nS^,  with  n reaching  62  at  t =1  and  n = 9 at  I = 99.  This  provides  complete  coverage  down  to 

100-sec  period.  We  have  toroidal  modes  Ttf,  with  I up  to  150  to  provide  complete  coverage 

ft  * 15 

down  to  60 -sec  period.  This  catalog  was  computed  by  Buland  and  Gilbert  and  was  kindly  made 
available  to  us. 


TABLE  111-6 

SOURCE  PARAMETERS  FOR  THE  ALASKAN  AND 

CHILEAN  EARTHQUAKES 

Parameter 

Alaska,  1964 

Chile,  1960 

Reference 

Depth  (km) 

50 

25 

16,  17 

Strike  (deg) 

246 

10 

16,  17 

Dip  (deg) 

20 

35.5 

16,  17 

Slip  (deg) 

270 

270 

16,  17 

Moment  (dyn-cm) 

7.5  X 1029 

8.6X  1029 

16,  17 

M 

s 

8.3 

8.3 

USGS;  USGS 

Energy  (erg) 

25 

1.5X  10  ° 

1.8X  1024 

16;  Gutenberg-Richter  relation 

As  a test  of  the  method,  we  consider  the  1964  Alaskan  and  I960  Chilean  earthquakes.  Both 

were  large  enough  to  excite  the  earth's  free  oscillations,  observations  of  which  were  used  to 

o 

determine  model  1066A.  The  source  mechanisms  of  these  were  found,  for  example,  by 
Kanamori1^  in  the  former  case,  and  by  Plafker  and  Savage17  in  the  latter.  The  pertinent  data 
are  listed  in  Table  III-6.  Double  couple  moment  tensors  computed  from  these  results  were 
then  used  along  with  the  catalog  of  normal  modes  to  calculate  the  energies  given  in  Table  III-7. 

The  results  for  both  earthquakes  show  that  the  step-function  source  releases  more  energy 
than  the  estimates  shown  in  Table  in-6;  actually,  it  can  be  shown  that  if  all  possible  modes  are 
considered,  the  energy  is  infinite  for  the  step  function.  Using  ramp  times  of  170  and  180  sec 
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TABLE  MI-7 

CALCULATED  ENERGIES  IN  THE  EARTH'S  FREE  OSCILLATIONS 

Toroidal 

Spheroidal 

Total 

from  Gutenberg-Richter 

(erg) 

(erg) 

(erg) 

Relation 

Alaska 

D = 45.2  km,  step  function 

2.55 X 1025 

25 

1.97X  10  3 

4.52  X 1025 

9.2 

D = 79.4  km,  step  function 

3.60  X 1025 

3.02  X 1025 

6.62  X 1025 

9.3 

D = 79.4  km,  170-sec  ramp  function 

3.39X  1023 

4.09  X 1024 

4.4  X1024 

8.5 

Chile 

D = 45.2  km,  180-sec  ramp  function 

3.35X  1023 

4.63X  1023 

7.98  X 1023 

8.1 

D = 45.2  km,  step  function 

2.64X  1025 

3.45  X 1025 

6.09  X 1025 

9.3 

for  the  Alaskan  and  Chilean  earthquakes,  respectively,  brings  both  energy  values  down  to  about 

a third  of  the  estimated  values.  These  ramp  times  were  suggested  by  analysis  of  other  large 

14 

earthquakes  and  seem  reasonable  when  compared  with  estimates  of  the  fault  length  divided  by 
16  18 

the  rupture  velocity.  ' Furthermore,  the  results  for  the  Alaskan  earthquake  with  the  source 
at  two  depths  show  a slight  increase  in  radiated  energy  with  increasing  depth. 

Two  conclusions  based  on  these  results  are  that:  (1)  although  the  seismic  moments  of  both 
the  Alaskan  and  Chilean  earthquakes  are  comparable,  the  difference  in  source  depth  and  mech- 
anism leads  to  five  times  more  energy  being  radiated  in  the  Alaskan  event;  and  (2)  the  bandwidth 
of  this  experiment  (100 -sec  period  for  the  spheroidal  modes  and  60 -sec  period  for  the  toroidal 
modes)  is  able  to  account  for  roughly  a third  of  the  estimated  strain  energy  release  with  a 170- to 
180-sec  ramp-function  time  history.  It  also  can  be  generally  observed  that  a very  large  part  of 
the  total  energy  is  contained  in  the  long  period  part  of  the  spectrum,  and  that  the  spectral  energy 
density  depends  significantly  on  the  source  time  function. 

Further  work  along  this  line  should  consist  of:  (1)  extending  the  spheroidal  part  of  the  free 
oscillation  catalog  to  shorter  periods,  (2)  analyzing  additional  events  to  determine  if  any  definite 
patterns  emerge  between  observed  source  characteristics  and  calculated  energies,  and  (3)  study- 
ing the  energy  distribution  as  a function  of  frequency. 

D.  W.  Me  Cow  an 
A.  M.  Dziewonski 

F.  APPLICATION  OF  NORMAL-MODE  THEORY  TO  CALCULATION  OF  CHANGES 
IN  THE  EARTH'S  MOMENT  OF  INERTIA  DUE  TO  EARTHQUAKES 

The  problem  of  evaluating  a change  in  the  earth's  moment  of  inertia  due  to  redistribution 
of  density  6p(r ) following  an  earthquake  consists  of  computation  of  the  integral: 

AC  = f [(r  • r)  I -(r  r)]  6p(r)  dV  . (HI-12) 

~ J y ~ 

In  all  published  solutions,  6p(  r ) has  been  obtained  by  solving  the  static  deformation  problem; 
here  we  use  the  result  of  Gilbert,  who  has  shown  that  the  static  displacement  field  can  be  ob- 
tained by  superposition  of  normal  modes.  In  practical  terms,  we  circumvent  the  difficulty  that 
occurs  in  solving  the  static  displacement  problem  with  the  "normal"  boundary  conditions  - for 

example,  continuity  of  the  radial  displacement  across  the  liquid-solid  interface  at  the  core 

19 

mantle  boundary.  Chinnery  makes  a convincing  case  for  application  of  a different  set  of  the 
boundary  conditions.  However,  in  the  case  of  a homogeneous,  adiabatically  compressed  fluid 
core  either  set  of  conditions  leads  to  the  same  system  of  equations.  Since  the  seismic  evidence 
indicates  that  the  density  distribution  in  the  outer  core  follows  the  Adams -William son  equation 
very  closely,  we  believe  that,  in  practice,  either  set  of  equations  should  lead  to  similar  results. 
In  particular,  in  our  numerical  calculations  we  use  an  earth  model^  that  satisfies  the  Adams- 
Williamson  equation  nearly  exactly. 

The  advantage  of  using  the  normal  mode  theory  is  that  it  has  been  empirically  tested;  all 
the  modes  that  significantly  contribute  to  the  change  in  the  moment  of  inertia  have  been  observed; 
their  excitation  and  eigenfrequencies  can  be  explained  by  the  realistic  models  of  the  earth's 
structure  and  seismic  source.  The  numerical  results  presented  by  various  authors  using  the 
static  displacement  approach  have  been  marred  by  disagreements  of  an  order  of  magnitude; 
comparisons  with  observations  are  questionable,  as  the  estimates  of  changes  in  the  pole  positions 
following  (or  preceding)  earthquakes  are  uncertain. 
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For  an  elastic  deformation: 


6p(  r)  = -pQ(r)  v • u(r ) -r  • u(r)arPc(r) 


(III  — 13) 


The  static  displacement  vector  u(r,t  = «>)  may  be  represented  in  terms  of  superposition  of  normal 
1 2 

modes  (the  epicentral  coordinate  system  is  adopted): 


u(£) = £ 

i=l 


«,  » min(2,i) 

ZEE 

1=0  n=0  m=0 


InC^i 


M. 

i 


(III- 14 ) 


where  M.  are  the  elements  of  the  moment  tensor  (the  formalism  and  symbols  follow  Gilbert  and 
1 9 

Dziewonski,  Sec.  2;  note  that  M.  is  a vector  representing  a second  rank  symmetric  tensor; 
similarly,  the  inertia  tensor  will  be  represented  as  a vector  C.). 

After  substitution  of  Eq.  (III-14)  into  (III-13)  and  then  into  (III-12),  transformed  into  the 
spherical  coordinate  system,  it  is  clear  that  each  element  of  inertia  tensor  must  be  a linear 
combination  of  the  elements  of  the  moment  tensor 


6 

iCi  - 2 v*j 

j=1 

where  D is  a matrix  representing  a fourth  rank  tensor: 


(III-15) 


and 


^1 


R2  - 2S2 


0 0 0 


R1  + S1  R2  + S2  S3  R2  + S2  + S3  0 0 


R1  + S1  R2  + S2  + S3  R2  + S2  S3  0 0 


s4  0 0 


0 s4  0 


0 0 — 2S, 


(III— 1 6) 


R1  = ! E nI0-  n(6i°)o 


n=0 


R2  = I E nI0  • n(62\ 


n=0 


1 = , nr  E nlz'  ; S2  “ 7"7?  E nh  ' n(S2^ 


3 \I~5 


n=0 


3 nT5 


n=0 


S3  = 75  E nh'n^zh  ’ S4  = </S  Z • <m-17> 


n=0 


n=0 
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The  are  defined  in  Eq.  (2.1.22)  of  Ref.  9 and 

J r3[2U(r)  + i(i  + 1)  V(r)]  po(r)  dr  . (III-18) 

Applicability  of  the  normal  mode  approach  depends  on  the  rate  of  convergence  of  the  infinite 
sums: 


i „i,  • n<«r>,  • (“-*  ’> 

n=0 

Since  the  strains  6 do  not  increase  with  n,  in  general,  the  convergence  depends  on  the 
behavior  of  the  integrals  nI  . Numerical  experiment  shows  that  the  sums  (in-19)  converge  very 
rapidly;  over  95  percent  of  the  entire  contribution  comes  from  only  three  modes:  0SQ,  QS2,  and 

1S2- 

Figure  III-4  shows  the  seven  distinct  elements  of  the  matrix  D as  a function  of  depth  for  a 
realistic  earth  model  PEM^  and  a uniform  earth  model^  with  p =5.52  g/cm^;  jl  = 1.463  Mb 
and  'X  = 3.514  Mb.  It  can  be  seen  that  some  functions,  such  as  — 2S^  and  R^  - 2S2,  have 
distinctly  different  behavior  for  the  two  earth  models.  It  follows  that  application  of  the  uniform 
earth  as  an  approximation  in  predicting  changes  in  the  moment  of  inertia  may  give  misleading 
results;  this  effect  will  depend  on  the  actual  source  mechanism. 

In  Table  III- 8 we  compare  results  obtained  for  several  source  mechanisms  with  the  values 
21 

of  Dahlen,  when  applicable.  The  agreement  is  excellent;  small  discrepancies  can  be  explained 
by  difference  in  the  earth  models  used  by  us  and  Dahlen. 

Most  of  the  large  earthquakes  occur  at  plate  boundaries;  it  has  been  established  before  that 
there  is  an  excellent  agreement  between  the  source  mechanism  inferred  from  excitation  of  seis- 
mic waves  and  that  predicted  by  the  relative  plate  motions. 

22 

Using  the  data  of  Solomon  et  al_.  for  the  positions  of  the  plate  boundaries  and  relative 
velocities  of  the  plates,  as  well  as  the  values  of  the  dip  of  Benioff  zones  determined  from  the 


TABLE  III— 8 

POLAR  SHIFT  (S;  in  units  of  0". 01 ) AND  ITS  DIRECTION  (D) 
FOR  SEVERAL  SOURCE  MECHANISMS  AND  EARTH  MODELS 

Earthquake 

Dahlen  (1973)21 

PEM-A 

PEM-C 

Uniform  Earth 

s 

D 

S 

D 

S 

D 

S 

D 

Chile17 

1.02 

1 10E 

1.39 

1 12E 

1.41 

1 12E 

0.24 

164E 

Chile18 

- 

- 

1.55 

1 19E 

1.56 

1 19E 

0.38 

178E 

Chile24 

- 

- 

3.45 

1 19E 

3.47 

1 19E 

0.84 

178E 

Alaska2^ 

0.48 

193E 

0.48 

194E 

0.48 

194E 

0.21 

156E 

Alaska1^ 

0.73 

202E 

0. 87 

204E 

0.88 

204 E 

0.38 

185E 
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earthquake  distribution,  we  were  able  to  estimate  fault  plane  parameters  for  hypothetical  earth- 
quakes of  three  types:  dip-slip  faulting  at  mid-oceanic  ridges,  strike-slip  faulting  on  fracture 
zones,  and  thrust  faulting  associated  with  the  subduction  zones.  In  Fig.  Ill- 5 we  show  the  magni- 
tude and  the  direction  of  the  change  in  the  pole  position  due  to  hypothetical  earthquakes  of  a 
standard  moment  of  10^  dyn-cm.  It  can  be  seen  that  there  is  a wide  range  of  the  potential  con- 
tribution of  different  seismic  zones  to  excitation  of  the  Chandler  wobble.  It  may  be  important 

that  the  arrows  associated  with  the  earthquakes  in  some  of  the  most  active  regions  of  the  circum- 

23 

Pacific  belt  have  approximately  the  same  direction.  Chinnery  and  Landers  present  evidence 
for  a simultaneous  (on  the  scale  of  1 to  2 months)  increase  of  seismicity  in  the  entire  Pacific 
region.  Thus,  contributions  of  earthquakes  in  different  zones  could  be  additive. 

A.  M.  Dziewonski 
R.  J.  O'Connellt 
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Fig.  Ill— 1 . Error  curve  (solid  curve) 
and  Akaike!s  FPE  (short  dashed  curve) 
and  AIC  (long  dashed  curve)  order  cri- 
terion curves  for  imaginary  part  of  log 
spectrum  indicating  that  data  have  an 
order  of  23. 
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Fig.  IE-2.  Raw  signals  at  subarray  centers  OC9  and  OC13  for  events  NTS  and  EK. 
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Fig.  IH-3.  Histograms  of  pP  detections  for  events  NTS  and  EK. 


Fig.  III-4.  Distinct  elements  of  matrix  Q 
[Eqs.  (Ill— 1 5)  and  (HI-16)]  as  a function  of 
depth  for  two  earth  models. 
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Fig.  Ill- 5.  Magnitude  (arrow's  length)  and  direction  (see  inset  for  definition) 
of  polar  shifts  in  response  to  three  types  of  hypothetical  earthquakes  with 
standard  moment  of  1030  dyn-cm. 


IV.  COMPUTER  SYSTEMS  AND  SOFTWARE 


A.  SRO  DATA  TAPE  PROCESSING  FACILITY 

PDP-7  Console  expansion  and  other  software  required  to  effectively  handle  SRO  data  tapes 
were  discussed  in  the  last  SATS.^  Several  of  the  design  and  implementation  tasks  mentioned 
have  now  been  completed. 

The  major  software  items  which  have  been  implemented  since  the  previous  report  are: 

(1)  Console  users  can  now  time- shift  and  modify  the  gain  of  an  individual 
seismic  trace  on  the  display  without  changing  any  of  the  other  traces. 

(2)  A PDP-11  program  (SROSP)  which  strips  short-period  data  from  an  SRO 
9 -track  tape  and  creates  a Data  Unit^  on  a 7 -track  tape  is  operational. 

Another  very  similar  program  to  do  the  same  for  the  long-period  SRO 
data  is  not  yet  implemented,  but  will  be  once  >DUI  and  >TAB  discussed 
below  have  been  completed. 

(3)  A PDP-7  program  (DUM)  to  merge  Data  Unit  tapes  generated  by  SROSP 
(and  the  not  yet  implemented  long-period  equivalent)  into  composite  Data 
Unit  tapes  containing  data  from  several  stations  is  now  available. 

(4)  A new  Console  program  (>SAV)  to  save  data  from  a Console  session  in 
the  same  Data  Unit  format  as  that  generated  by  SROSP  and  DUM  has 
been  designed  and  completed.  It  will  append  data  to  an  existing  Data 
Unit  as  well  as  create  a new  one. 

(5)  As  part  of  the  above  software  implementation,  several  general  Fortran 
subroutines  were  created  for  our  PDP-7  computers.  A set  of  these 
allow  character  string  manipulations.  Others  create,  check,  and  ex- 
tract information  from  the  Type,  Description,  and  Title  sections  of  a 
Data  Unit.  Still  others  support  reading  and  writing  of  the  actual  data 
in  the  body  of  the  Data  Unit. 

(6)  Also  in  support  of  this  software  effort,  the  PDP-7  Fortran  system  was 
modified  to  allow  multiple  formatted  reads  of  the  same  physical  record. 

Other  software  components  have  been  designed  and  are  now  being  coded,  the  major  one  of 
which  is  >DUI  which  will  provide  convenient  and  efficient  initialization  of  the  Console  from  Data 
Unit  tapes.  This  program  will  be  driven  by  a condition  table  set  up  by  another  program,  >TAB. 
The  first  version  of  >TAB  is  also  specified.  It  will  set  up  the  condition  table  to  find  all  SRO 
P-  or  Rayleigh-wave  data  for  an  event,  given  only  the  event  location  and  origin  time. 

It  will  be  inconvenient  for  researchers  to  handle  multistation  SRO  data  in  the  Console  envi- 
ronment until  >DUI  and  >TAB  are  completed.  However,  a simple  test  program  is  now  available 
which  does  allow  us  to  enter  Data  Units  into  the  Console.  For  example.  Fig.  IV- 1 shows  twenty 
SRO  short-period  waveforms  which  were  entered  into  the  Console  and  then  plotted  with  existing 
programs.  Table  IV- 1 lists  events  reported  by  USGS  which  correspond  to  some  of  the  signals 
in  Fig.  IV- 1. 


t Seismic  Discrimination  SATS,  Lincoln  Laboratory,  M.I.T.  (30  June  1975),  DDC  AD-A014793/4. 
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TABLE  IV-1 

ASSOCIATION  OF  SRO  DATA  TRACES  OF  FIGURE  IV-1 
WITH  EVENTS  REPORTED  BY  USGS 

Event  Parameters 

Trace 

Number 

Date 

(1974) 

Origin 

Time 

Latitude 

(°S) 

Longitude 

(°W) 

Magnitude 

Distance 

(deg) 

1 

22  December 

02:18:36.9 

26.9 

176.2 

5.1 

90.1 

8(1) 

22  December 

16:44:15.5 

19.0 

69.8 

5.2 

64.0 

8(2) 

22  December 

16:44:05.3 

17.6 

179.0 

5.1 

86.3 

10 

22  December 

18:06:36.2 

27.0 

176.3 

5.2 

90.3 

15 

23  December 

01:04:02.7 

14.6 

175.7 

5.4 

82.0 

17 

23  December 

08:13:58.8 

3.0 

75.4 

4.6 

47.7 

18 

23  December 

11:16:48.2 

16.2 

176.7 

4.6 

83.8 

19 

23  December 

11:31:00.0 

26.3 

177.5 

4.2 

90.7 

(1),(2)  Arrivals  from  two  events  are  on  the  same  trace.  First  arrivals  for  each  are  indicated  by  (1) 
and  (2)  on  trace  8 of  Fig.  IV-1. 

Each  SRO  contains  a short-period  event  detector  which  triggers  the  saving  of  the  data.  The 
twenty  waveforms  shown  in  Fig.  IV-1  correspond  to  the  first  twenty  detections  by  the  AIjQ  site 
starting  on  22  December  1974.  All  gains  are  the  same  except  traces  15  and  20  which  are  atten- 
uated by  1/4  and  1/8,  respectively.  Those  signals  which  correspond  to  events  on  the  USGS 
Preliminary  Determination  of  Epicenters  List  are  identified  in  Table  IV-1.  In  one  case  (trace  8) 
two  events  occur,  with  the  second  one  arriving  about  65  sec  after  the  first.  Traces  5,  9,  and 
16  show  no  clear  arrivals  but  may.  contain  coda  or  secondary  arrivals  from  the  events  with  the 
first  arrivals  on  traces  4,  8,  and  15,  respectively.  There  were  four  PDE  events  in  the  same 
time  period  with  distances  <90°  from  ALQ  which  did  not  trigger  the  ALQ  SRO  event  detector. 

The  PDE  magnitudes  and  distances  for  these  events  are: 


^b 

Distance 

(deg) 

5.0 

70 

4.6 

53 

3.5 

69 



39 

Several  of  the  traces  not  associated  with  the  PDE  list  show  very  clear  arrivals  and,  as 
mentioned  above,  some  do  not  strongly  indicate  distinct  arrivals.  Trace  3,  however,  shows  a 
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relatively  clear  but  small  arrival.  Figure  IV-2  shows  trace  3 and  a filtered  version  which 
clearly  shows  the  arrival.  The  filtered  version  was  obtained  simply  by  the  application  of  the 
Console  filter  program. 

The  SBO  at  ALQ  was  just  starting  operation  in  December  1974.  It  is  not  our  purpose  here 
to  evaluate  the  site  or  the  performance  of  the  event  detector  at  that  time.  However,  some  ob- 
servations can  be  made  about  our  small  data  set.  The  data  quality  is  excellent.  The  event  de- 
tector was  set  in  such  a way  that  it  generated  very  few  false  alarms.  It  might  be  that  if  more 
false  alarms  had  been  allowed,  one  or  two  of  the  missed  PDE  events  would  have  triggered  the 
detector.  Certainly,  in  a network  environment  many  more  false  alarms  can  be  and  should  be 
tolerated  than  occur  in  our  small  sample  of  data.  ^ ^ Lacoss 

M.  F.  O’Brien 
L.  J.  Turek 

B.  USE  OF  ARPANET  FOR  INTERACTIVE  DISPLAY  AND  ANALYSIS 

OF  SEISMIC  NETWORK  DATA 

A study  of  earthquake  source  mechanisms  using  P-waves  from  shallow  Asian  earthquakes 
was  reported  in  the  last  SATS  (30  June  1975,  DDC  AD-A014793/4).  That  study  involved  the  use 
of  seismograms  from  a network  of  stations  and  of  theoretical  seismograms  for  the  same  sta- 
tions. Displays  showing  theoretical  and  observed  seismograms  and  how  the  station  was  located 
with  respect  to  the  hypocenter  and  fault  plane  orientations  were  found  to  be  quite  helpful  in  that 
study.  However,  the  construction  of  the  displays  was  a slow  manual  process  which  detracted 
from  the  study  and  use  of  the  data.  Motivated  by  this,  an  interactive  seismic  display  and  analy- 
sis package  has  been  written  to  generate  displays  and  perform  the  kind  of  analysis  required  for 
those  source-mechanism  studies.  The  package  is  a fairly  general  interactive  display  system 
which  has  been  implemented  to  demonstrate  how  remote  computation  and  data  storage  can  be 
combined  with  a graphics  terminal  as  a tool  for  seismic  research  involving  a number  of  seismic 
stations.  The  source-mechanism  problem  is  just  one  possible  application. 

The  display  and  analysis  package  runs  under  the  time-sharing  system  on  the  Lincoln  Lab- 
oratory IBM-370/168.  The  user  terminal  is  a Tektronics  terminal  with  a storage  cathode-ray- 
tube  display.  The  terminal  is  attached  to  the  370  through  the  ARPANET.  Two  such  terminals 
are  available  and  are  used  for  a variety  of  interactive  and  graphics-based  activities. 

Figure  IV- 3 is  an  example  of  the  kind  of  display  generated  by  the  system.  It  can  display  up 
to  twenty  waveforms  for  an  event,  the  orientation  of  the  stations  with  respect  to  the  event  loca- 
tion, and  the  orientation  with  respect  to  the  fault  plane  and  slip  vector  if  these  are  known.  In 
the  figure  shown,  the  central  circle  is  the  lower  half  focal  sphere  for  the  event.  A symbol  is 
shown  in  that  circle  for  each  station  of  interest.  The  symbol  is  at  the  azimuth  of  the  station  from 
the  event,  and  the  distance  from  the  origin  represents  the  takeoff  angle  for  P-waves.  The  fault 
plane  is  also  shown  in  the  usual  way  on  such  a display.  The  symbol  used  for  each  station  indi- 
cates if  the  first  motion  is  up,  down,  or  has  not  been  entered  into  the  data.  More  stations  can 
be  displayed  than  seismograms  and,  in  fact,  the  display  can  be  used  for  fault-plane  studies 
without  displaying  any  seismograms  at  all.  In  the  display  shown,  there  are  two  seismograms 
shown  for  each  station.  One  of  these  is  the  actual  seismogram  of  the  long-period  P-wave,  and 
the  other  is  a theoretical  one  generated  by  the  model  built  into  the  system.  In  each  case,  the 
actual  seismogram  is  identified  by  the  station  three-letter  code  and  the  theoretical  one  has  a 
T appended.  The  * indicates  that  time  scaling  has  been  adjusted  and  the  user  should  refer  to 
the  display  system  output  file,  which  is  discussed  below,  for  details. 
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The  display  and  analysis  package  makes  use  of  four  distance  files  of  data  - three  are  input 
and  one  is  output.  Two  of  the  input  files  constitute  a waveform  library  which  can  be  created 
and  modified  independent  of  the  display  system.  These  files  can  contain  any  number  of  seismo- 
grams. They  can  be  both  long  and  short  period,  actual,  or  theoretical.  Each  is  identified  by 
station,  date,  time,  time  interval,  sampling  rate,  gain,  instrument  type,  a comment  if  desired, 
and  a unique  seismogram  identification  number.  These  two  files  are  the  universe  of  possible 
signals  which  can  be  displayed  by  the  display  system.  If  the  display  system  calls  a subsystem 
to  generate  theoretical  waveforms,  the  theoretical  ones  are  added  to  this  universe  and  treated 
exactly  like  actual  signals  for  display  purposes.  The  third  input  file  is  specific  to  the  display 
package.  It  contains  the  hypocenter  and  origin  time  of  events  of  interest  and,  if  desired,  speci- 
fications of  the  stations  and  waveforms  of  interest.  If  such  specifications  are  not  given,  the 
program  will  immediately  enter  the  interactive  mode  and  display  a menu  of  options  to  the  user. 
The  output  file  just  contains  tables  which  completely  specify  the  actual  display.  In  fact,  the 
display  is  generated  by  internal  subroutines  which  are  driven  by  this  so-called  output  file  and 
which  get  waveform  data  from  the  so-called  input  files. 

The  display  shown  in  Fig.  IV- 3 was,  in  fact,  generated  using  the  interactive  capability  of 
the  system.  The  original  display  generated  for  the  same  event  is  shown  in  Fig.  IV-4.  During 
the  interactive  session,  the  user  has  generated  theoretical  seismograms  and  added  them  to  the 
display,  changed  various  scaling  parameters,  added  the  fault  plane,  and  changed  the  size  of 
the  central  circle  to  clean  up  the  general  appearance  of  the  plot. 

Table  IV-2  shows  the  menu  available  to  the  user.  At  essentially  any  point  in  a terminal 
session,  he  can  request  to  exercise  any  or  all  of  the  options.  He  simply  enters  the  list  of 


MENU  OF 

TABLE  IV-2 

OPTIONS  FOR  THE  INTERACTIVE  DISPLAY  SYSTEM 

Argument 

Definition 

WAVFM 

Display  waveform  file  before  making  a change 

ISEIS 

Add,  replace,  or  delete  waveforms  to  display 

SCLRAD 

Scale  the  radius  of  the  net  projection  to  another 

MAXPTS 

Change  maximum  number  of  waveform  points  to  display 

FAULT 

Display  the  fault  plane  and  theoretical  waveforms 

CNET 

Change  the  stereographic  projection  to  azimuth-delta  type 

CTOA 

Change  the  delta  limits  for  the  station  projections 

CONLIN 

Do  not  connect  plot  projection  and  waveform  with  a line 

NOAZI 

Do  not  order  waveforms  azimuthally  for  plot 

CTIM 

Scale  the  time  axis  for  the  waveform(s) 

CAMP 

Scale  the  amplitude  display  for  the  waveform(s) 

CSTA 

Change  the  station  code(s)  to  waveform  serial  number(s) 

TMSTRT 

Change  the  start  time  of  the  displayed  waveform(s) 
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arguments  corresponding  to  the  changes  he  wishes  to  make,  and  the  system  then  requests  all 
further  information  it  requires.  The  normal  use  of  the  system  is  to  set  an  initial  display,  make 
a hard  copy  using  a device  attached  to  the  terminal,  request  the  menu  and  make  changes, 
generate  a new  display,  etc.  L C I ande 

R.  T.  Lacoss 


C.  PROGRESS  REPORT  ON  NLS  (ON-LINE  SYSTEM) 

AND  THE  SEISMIC  DATA  MANAGEMENT  SYSTEM 

Most  of  the  major  files  that  will  make  up  the  NLS  directory  for  the  Seismic  Data  Manage- 
ment System  have  been  updated  and  structured  within  the  NLS  system.  Some  of  these  files  con- 
tain the  site  descriptions  of  the  High  Gain  Long  Period  sites  (HGLP)  and  of  the  Seismic  Research 
Observatories  (SROs).  They  reflect  the  current  state  of  development  in  those  areas  where  de- 
velopment is  still  in  progress,  the  Seismic  Research  Observatories  for  example,  as  well  as  con- 
taining complete  descriptions  of  both  the  hardware  aspect  of  the  SROs  and  the  software  necessary 
to  run  the  sites. 

The  main  effort  that  is  now  under  way  is  basically  divided  into  two  parts.  The  first  part  is 
the  continuing  acquisition  of  additional  information  and  documentation  that  will  be  placed  in  the 
directory  <IWWSS>.  Information  related  to  the  Network  Event  Processor  (NEP),  Seismic  In- 
formation Processor  (SIP),  Control  and  Communication  Processor  (CCP),  as  well  as  the  actual 
structure  of  the  seismic  files  that  will  be  stored  on  the  mass  store  device,  would  fall  into  the 
category  of  information  that  should  be  contained  in  the  NLS  directory.  Also  connected  with  this 
aspect  of  the  effort  is  the  continued  acquisition  of  updated  information  describing  the  installa- 
tion of  the  SRO  sites,  and  the  additional  information  about  the  SRO  installations  as  they  become 
operational  and  begin  to  produce  actual  data. 

The  second  main  area  of  effort  concerns  procedures  by  which  researchers  and  other  users 
of  the  seismic  information  directory  will  access  this  information.  All  the  information  that  is 
contained  and  collected  in  the  seismic  directory  <IWWSS>  will  be  made  available  in  hard-copy 
form,  and  therefore  not  require  that  the  user  have  any  prior  knowledge  of  the  ARPANET  or  NLS. 
This  is  the  simplest  method  of  getting  the  information  to  interested  users.  However,  this  is  a 
slow,  cumbersome,  and  inappropriate  means  of  information  transfer.  The  other  means  avail- 
able are  to  provide  access  to  the  directory  <IWWSS>  and  allow  the  user  to  obtain  those  files  that 
pertain  to  the  topics  of  interest.  The  files  maintained  in  this  directory  will  be  in  both  text  form 
and  in  NLS  form.  The  advantage  to  the  user  is  that  he  should  be  able  to  copy  these  files  on  his 
terminal  as  a sequential  text  file,  or  he  can  enter  the  NLS  world  and  freely  browse  through  these 
NLS  files  addressing  any  topics  of  interest. 

There  are  two  approaches  to  providing  on-line  ARPANET  access  to  the  data  files.  The 
simplest  and  most  convenient  approach  would  be  to  make  these  files  available  in  a semipublic, 
protected  directory  with  read-only  privilege  similar  to  the  NIC  directory  at  OFFICE- 1.  This 
is  the  approach  that  is  now  being  taken.  The  directory  <IWWSS>  will  have  the  proper  protection 
and  restrictions  to  allow  the  user  access  to  only  this  directory  and  not  the  rest  of  the  OFFICE-1 
(or  ISIC)  system.  The  data  information  files  would  be  maintained  both  as  NLS  files  and  as  text 
files,  thus  allowing  the  user  the  choice  of  methods  of  data  retrieval.  The  text  versions  could  be 
listed  at  the  terminal  of  the  user  or  at  some  other  line  printer,  while  the  NLS  versions  would 
be  available  to  experienced  NLS  users. 
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The  developers  of  NLS,  the  Stanford  Research  Institute,  have  now  begun  to  train  and  dis- 
tribute information  which  will  allow  users  to  write  and  set  up  complete  user  subsystems  within 
the  NLS  environment.  This  will  provide  a very  useful  and  powerful  means  for  setting  up  a com- 
plete seismic  user  subsystem  (SUS)  within  NLS.  The  complete  design  of  SUS  is  not  yet  complete; 
however,  the  goal  of  this  subsystem  will  be  to  interact  with  the  user  and  try  to  provide  him  with 
description  and  information  that  is  desired  without  making  a novice  user  stumble  through  the 
data  files.  R.  M<  sheppard 


D.  THE  APPLIED  SEISMOLOGY  GROUP  (ASG) 

ARPANET  CONNECTION 

The  ASG  uses  the  ARPANET  to  provide  time-sharing  terminal  access  to  several  large  com- 
puter sites,  as  described  in  the  previous  SATS  (30  June  1975,  DDC  AD-A014793/4).  An  ELF 
operating  system,  running  on  a PDP- 11/40  computer,  connects  seven  terminals  and  a lineprinter 
to  the  ARPANET. 

We  are  continuing  to  improve  the  software  and  the  hardware  supporting  the  network 
connection.  Recent  improvements  include: 

Improved  reliability  The  ELF  system  crashes  two  times  per  week,  on  the 

average.  This  represents  an  order-of-magnitude  im- 
provement in  reliability  from  that  previously  reported. 

The  ELF  system  was  transferred  to  the  PDP-ll/40, 
releasing  the  11/50  (which  has  many  peripheral  de- 
vices) for  software  and  hardware  development. 

We  have  designed  and  tested  a simple  but  robust  high- 
speed interprocessor  device  for  PDP-11  computers. 

Three  pairs  of  these  will  be  built,  enabling  us  to 
communicate  efficiently  among  our  three  PDP-ll's 
(an  11/50,  an  11/40  and  a GT-40  graphics  terminal). 

Software  provided  on  the  Lincoln  370  computer  by 
Group  28  to  allow  sending  files  over  the  ARPANET 
can  now  be  used  to  send  files  to  our  ELF  system  for 
lineprinter  listing.  Previously,  we  had  been  able  to 
get  listings  only  from  TENEX  sites. 

P.  A.  Neilson 
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Fig.  IV-2.  Unfiltered  short-period  vertical  SRO  waveform,  and  same  trace 
bandpass- filtered  to  pass  0.6  to  3.0  Hz. 
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Fig.  IV- 3.  Actual  and  theoretical  seismograms  displayed 
by  interactive  display  system. 
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Fig.  IV- 4.  Initial  seismogram  display  from  interactive 
display  system. 
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